This document can be cited as: Draine, B.T., and Flatau, PJ. 2009, 
"User Guide for the Discrete Dipole Approximation Code DDSCAT 7.0", 
http:7/arxiv.o rg/abs/0809. 0337v5 1 



User Guide for the Discrete Dipole 
Approximation Code DDSCAT 7.0 

Bruce T. Draine 
Princeton University Observatory 
Princeton NJ 08544-1001 

(draine@ astro .princeton. edu) 
and 

Piotr J. Flatau 
University of California San Diego 
Scripps Institution of Oceanography 
La JoUa CA 92093-0221 

(pf latau@ucsd . edu) 

last revised: 2009 March 9 



Abstract 

DDSCAT 7.0 is a freely available open-source Fortran-90 software package applying the "discrete 
dipole approximation" (DDA) to calculate scattering and absorption of electromagnetic waves by targets 
with arbitrary geometries and complex refractive index. The targets may be isolated entities (e.g., dust 
particles), but may also be 1-d or 2-d periodic arrays of "target unit cells", which can be used to study 
absorption, scattering, and electric fields around arrays of nanostructures. 

The DDA approximates the target by an array of polarizable points. The theory of the DDA and its 
implementation in DDSCAT is presented in Draine (1988) and Draine & Flatau (1994), and its extension 
to periodic structures (and near-field calculations) in Draine & Flatau (2008). DDSCAT 7.0 allows 
accurate calculations of electromagnetic scattering from targets with "size parameters" 27racfr/A ^ 25 
provided the refractive index m is not large compared to unity (|m — 1| ^E, 2). DDSCAT 7.0 includes 
support for MP I, OpenMP, and the Intel® Math Kernel Library (MKL). 

DDSCAT supports calculations for a variety of target geometries (e.g., ellipsoids, regular tetrahedra, 
rectangular solids, finite cylinders, hexagonal prisms, etc.). Target materials may be both inhomoge- 
neous and anisotropic. It is straightforward for the user to "import" arbitrary target geometries into the 
code. DDSCAT automatically calculates total cross sections for absorption and scattering and selected 
elements of the Mueller scattering intensity matrix for specified orientation of the target relative to the 
incident wave, and for specified scattering directions. DDSCAT 7.0 can calculate scattering and absorp- 
tion by targets that are periodic in one or two dimensions. DDSCAT 7.0 also supports postprocessing to 
calculate E and B at user-selected locations near the target; a Fortran-90 code DDfleld for this purpose 
is included in the distribution. 

This User Guide explains how to use DDSCAT 7.0 to carry out electromagnetic scattering calcula- 
tions. 
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1 Introduction 

DDSCAT is a software package to calculate scattering and absorption of electromagnetic waves by tar- 
gets with arbitrary geometries using the "discrete dipole approximation" (DDA). In this approximation 
the target is replaced by an array of point dipoles (or, more precisely, polarizable points); the electromag- 
netic scattering problem for an incident periodic wave interacting with this array of point dipoles is then 
solved essentially exactly. The DDA (sometimes referred to as the "coupled dipole approximation") was 
apparently first proposed by Purcell & Pennypacker (1973). DDA theory was reviewed and developed 
further by Draine (1988), Draine & Goodman (1993), reviewed by Draine & Flatau (1994) and Draine 
(2000), and recently extended to periodic structures by Draine & Flatau (2008). 

DDSCAT 7.0, the current release of DDSDCAT, is an open-source Fortran 90 implementation of the 
DDA developed by the authorsQ It extends DDSCAT to treat structures that are periodic in one or two 
dimensions, using methods described by Draine & Flatau (2008). 

DDSCAT is intended to be a versatile tool, suitable for a wide variety of applications including 
studies of interstellar dust, atmospheric aerosols, blood cells, marine microorganisms, and nanostructure 
arrays. As provided, DDSCAT 7.0 should be usable for many applications without modification, but the 
program is written in a modular form, so that modifications, if required, should be fairly straightforward. 

The authors make this code openly available to others, in the hope that it will prove a useful tool. We 
ask only that: 

• If you publish results obtained using DDSCAT, please acknowledge the source of the code, and 
cite relevant papers, such as Draine & Flatau (1994), Draine & Flatau (2008), and this UserGuide 
(Draine & Flatau 2008). 

• If you discover any errors in the code or documentation, please promptly communicate them to the 
authors. 

• You comply with the "copyleft" agreement (more formally, the GNU General Public License) of 
the Free Software Foundation: you may copy, distribute, and/or modify the software identified as 
coming under this agreement. If you distribute copies of this software, you must give the recipients 
all the rights which you have. See the file doc /copyleft . txt distributed with the DDSCAT 
software. 

We also strongly encourage you to send email to draine@astro.princeton.edu identifying 
yourself as a user of DDSCAT; this will enable the authors to notify you of any bugs, corrections, or 
improvements in DDSCAT. Up-to-date information on DDSCAT can be found at 

| http : / / www . astro _. princeton . edu/ ~dr a ine/ DDSCAT . html| 
Information is also available at http : / /ddsc at . wikido t . com/start| 

The current version, DDSCAT 7.0, uses the DDA formulae from Draine (1988), with dipole polariz- 
abiUties determined from the Lattice Dispersion Relation (Draine & Goodman 1993, Gutkowicz-Krusin 
& Draine 2004). The code incorporates Fast Fourier Transform (FFT) methods (Goodman, Draine, & 
Flatau 1991). DDSCAT 7.0 includes capability to calculate scattering and absorption by targets that are 
periodic in one or two dimensions - arrays of nanostructures, for example. The theoretical basis for 
application of the DDA to periodic structures is developed in Draine & Flatau (2008). 

' The release history of DDSCAT is as follows: 

• DDSCAT 4b: Released 1993 March 12 

• DDSCAT 4bl: Released 1993 July 9 

• DDSCAT 4c: Although never announced, DDSCAT . 4 c was made available to a number of interested users beginning 1994 
December 18 

• DDSCAT 5a7: Released 1996 

• DDSCAT 5a8: Released 1997 April 24 

• DDSCAT 5a9: Released 1998 December 15 

• DDSCAT 5al0: Released 2000 June 15 

• DDSCAT 6.0: Released 2003 September 2 

• DDSCAT 6.1: Released 2004 September 10 

• DDSCAT 7.0: Released 2008 September 1 
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We refer you to the list of references at the end of this document for discussions of the theory and 
accuracy of the DDA [in particular, reviews by Draine and Flatau (1994) and Draine (2000), and recent 
extension to 1-d and 2-d arrays by Draine & Flatau (2008).]. In Sj4]we describe the principal changes 
between DDSCAT 7.0 and the previous releases. The succeeding sections contain instructions for: 

• compiling and Unking the code; 

• running a sample calculation; 

• understanding the output from the sample calculation; 

• modifying the parameter file to do your desired calculations; 

• specifying target orientation; 

• changing the DIMENSIONing of the source code to accommodate your desired calculations. 

The instructions for compiling, linking, and running will be appropriate for a Linux system; slight 
changes will be necessary for non-Linux sites, but they are quite minor and should present no difficulty. 

Finally, the current version of this User Guide can be obtained from 
|http ://arxiv.org/abs/0809.0337| - you will be offered the options of downloading either 
Postscript or PDF versions of the User Guide. 



2 Applicability of the DDA 

The principal advantage of the DDA is that it is completely flexible regarding the geometry of the target, 
being limited only by the need to use an interdipole separation d small compared to (1) any structural 
lengths in the target, and (2) the wavelength A. Numerical studies (Draine & Goodman 1993; Draine & 
Flatau 1994; Draine 2000) indicate that the second criterion is adequately satisfied if 

\m\kd<l , (1) 

where m is the complex refractive index of the target material, and k = 2tt/X, where A is the wavelength 
invacuo. This criterion is valid provide |7ti— 1| ^ 3orso. When Im(m) becomes large, the DDA solution 
tends to overestimate the absorption cross section Cabs, and it may be necessary to use interdipole 
separations d smaller than indicated by eq. ([TJ to reduce the errors in Cabs to acceptable values. 

If accurate calculations of the scattering phase function (e.g., radar or lidar cross sections) are desired, 
a more conservative criterion 

|m|fcd<0.5 (2) 

will usually ensure that differential scattering cross sections dCscn/dft are accurate to within a few 
percent of the average differential scattering cross section Csra/47r (see Draine 2000). 

Let V be the actual volume of solid material in the target|3 If the target is represented by an array of 
N dipoles, located on a cubic lattice with lattice spacing d, then 

V = Nd^ . (3) 

We characterize the size of the target by the "effective radius" 

fleff = (3F/47r)i/3 , (4) 

the radius of an equal volume sphere. A given scattering problem is then characterized by the dimen- 
sionless "size parameter" 

- , 27racff - 

X = /Cflcff = — r — . (5) 



A 



The size parameter can be related to N and \m\kd: 



27raeff 62.04 / iV ^ 



A ItoI V 10*^ 



• \m\kd . (6) 



^ In the case of an infinite periodic target, V is the volue of solid material in one "Target Unit Cell", 
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Equivalently, the target size can be written 



106 



Oeff = 9.873— — • \m\kd . (7) 



Practical considerations of CPU speed and computer memory currently available on scientific worksta- 
tions typically limit the number of dipoles employed to < 10^ (see fJT6]for limitations on due to 
available RAM); for a given N, the limitations on |m|fc(i translate into limitations on the ratio of target 
size to wavelength. 

For calculations of total cross sections Cabs and Csca, we require |m|fcd < 1: 

A / N 62.04 / N ^^^^ 



For scattering phase function calculations, we require \m\kd < 0.5: 

A / N 31.02 / N ^^^^ 



It is therefore clear that the DDA is not suitable for very large values of the size parameter x, or very 
large values of the refractive index m. The primary utility of the DDA is for scattering by dielectric 
targets with sizes comparable to the wavelength. As discussed by Draine & Goodman (1993), Draine 
& Flatau (1994), and Draine (2000), total cross sections calculated with the DDA are accurate to a few 
percent provided N > 10"' dipoles are used, criterion ([T]) is satisfied, and the refractive index is not too 
large. 

For fixed Imj/cd, the accuracy of the approximation degrades with increasing |m — 1|, for reasons 
having to do with the surface polarization of the target as discussed by CoUinge & Draine (2004). With 
the present code, good accuracy can be achieved for |to — 1 1 < 2. 

Examples illustrating the accuracy of the DDA are shown in Figs.[T]-(2]which show overall scattering 
and absorption efficiencies as a function of wavelength for different discrete dipole approximations to a 
sphere, with N ranging from 304 to 59728. The DDA calculations assumed radiation incident along the 
(1,1,1) direction in the "target frame". Figs.[3}Slshow the scattering properties calculated with the DDA 
for X ~ ka = 7. Additional examples can be found in Draine & Flatau (1994) and Draine (2000). 



3 DDSCAT 7.0 

3.1 What Does It Calculate? 

3.1.1 Absorption and Scattering by Finite Targets 

DDSCAT 7.0, like previous versions of DDSCAT, solves the problem of scattering and absorption by 
a finite target, represented by an array of polarizable point dipoles, interacting with a monochromatic 
plane wave incident from infinity. DDSCAT 7.0 has the capability of automatically generating dipole 
array representations for a variety of target geometries (see can also accept dipole array repre- 

sentations of targets supplied by the user (although the dipoles must be located on a cubic lattice). The 
incident plane wave can have arbitrary elliptical polarization (see ^2% . and the target can be arbitrarily 
oriented relative to the incident radiation (see 23- The following quantities are calculated by DDSCAT 
7.0: 

• Absorption efficiency factor Qabs = C'abs/'i'acff' where Cabs is the absorption cross section; 

• Scattering efficiency factor Qsca = Csca/Tra^jf , where Csca is the scattering cross section; 

• Extinction efficiency factor Qoxt = Qsca + Qahs', 

• Phase lag efficiency factor Qpha, defined so that the phase-lag (in radians) of a plane wave after 
propagating a distance L is just ntQphaT^alffL, where is the number density of tai^gets. 
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1 2 3 4 5 6 7 8 9 10 11 12 13 

x = ka 



Figure 1: Scattering and absorption for a spliere with m — 1.33 + O.Oli. Tlie upper panel shows the exact values of Qsca and 
Qabs, obtained with Mie theory, as functions of a; = ka. The middle and lower panels show fractional errors in Qsca and Qabs, 
obtained using DDSCAT with polarizabilities obtained from the Lattice Dispersion Relation, and labelled by the number A'^ of 
dipoles in each pseudosphere. After Fig. 1 of Draine & Flatau (1994). 
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x = ka 



Figure 2: Same as Fig.[Il but for ?n = 2 + z. After Fig. 2 of Draine & Flatau (1994). 
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Figure 3: Differential scattering cross section for m — 1.33 + O.Oli pseudosphere and ka = 7. Lower panel shows fractional 
error compared to exact Mie theory result. The A'^ — 17904 pseudosphere has \m\kd — 0.57, and an rms fractional error in da/dQ 
of 2.4%. After Fig. 5 of Draine & Flatau (1994). 




Figure 4: Same as Fig.[3]but for m = 2 + i. The A'^ = 59728 pseudosphere has \m\kd = 0.65, and an rms fractional error in 
da/dD. of 6.7%. After Fig. 8 of Draine & Flatau (1994). 
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• The 4x4 Mueller scattering intensity matrix Sij describing the complete scattering properties of 
the target for scattering directions specified by the user (see ^2M . 

• Radiation force efficiency vector Qiad (see i^TSl l. 

• Radiation torque efficiency vector Qr (see ifTS] ). 

In addition, the user can choose to have DDSCAT 7.0 store the solution for post-processing. A separate 
Fortran-90 program DDf ield (see f|27|i is provided that can readily calculate E and B at user-selected 
locations. 

3.1.2 Absorption and Scattering by Periodic Arrays of Finite Structures 

DDSCAT 7.0 includes the capability to solve the problem of scattering and absorption by an infinite 
target consisting of a 1-d or 2-d periodic array of finite structures, illuminated by an incident plane wave. 
The finite structures are themselves represented by arrays of point dipoles. 

The electromagnetic scattering problem is formulated by Draine & Flatau (2008), who show how the 
problem can be reduced to a finite system of linear equations, and solved by the same methods used for 
scattering by finite targets. 

The far-field scattering properties of the 1-d and 2-d periodic arrays can be conveniently represented 
by a generalization of the Mueller scattering matrix to the 1-d or 2-d periodic geometry - see Draine & 
Flatau (2008) for definition of S'fj'^^ (M, C) and S^^"^^ {M, N). 

As for finite targets, the user can choose to have DDSCAT 7.0 store the calculated polarization field 
for post-processing. A separate Fortran-90 program DD field is provided to calculate E and B at 
user-selected locations, which can be located anywhere (including within the target). 

3.2 Application to Targets in Dielectric Media 

Let oj be the angular frequency of the incident radiation. For many applications, the target is essen- 
tially in vacuo, in which case the dielectric function e which the user should supply to DDSCAT is the 
actual complex dielectric function etargot {^)^ or complex refractive index Wtargct (^) = Retarget of the 
target material. 

However, for many applications of interest (e.g., marine optics, or biological optics) the "target" body 
is embedded in a (nonabsorbing) dielectric medium, with (real) dielectric function emcdium('jj), or (real) 
refractive index m,„cdium('^) = -\/emodium- DDSCAT is fully applicable to these scattering problems, 
except that: 

• The "dielectric function" or "refractive index" supplied to DDSCAT should be the relative dielec- 
tric function 

or relative refractive index: 

m(w) 

• The wavelength A specified in ddscat . par should be the wavelength in the medium: 



^^mcdium 

where Avac = 2iic/ijj is the wavelength in vacuo. 

Because DDSCAT will be using the relative refractive index m{uj) and the wavelength in the medium, 
the file supplying the refractive index as a function of wavelength should give the relative refractive 
index or relative dielectric function as a function of wavelength in the medium. Those using DDSCAT 
for targets that are not in vacuo must be attentive to this when preparing the input table(s) giving the 
relative dielectric function or refractive index of the target material(s). 

The absorption, scattering, extinction, and phase lag efficiency factors Qabs. Qsca. and Qc^t cal- 
culated by DDSCAT will then be equal to the physical cross sections for absorption, scattering, and 



etargct(w) ^jQ-^ 



'^medium V 



"Itargot(w) 



m-modium(t^) 



(11) 
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extinction divided by Tra^fj . For example, the attenuation coefficient for radiation propagating through a 
medium with a number density rtf of scatterers will be just a — ritQcxtTra^g . Similarly, the phase lag 
(in radians) after propagating a distance L will be utQphn'^alffL. 

The elements Sij of the 4x4 Mueller scattering matrix S calculated by DDSCAT for finite targets 
will be correct for scattering in the medium: 



-) 

27rr J 



2 

S • lin, (13) 



where Ijn and Igca are the Stokes vectors for the incident and scattered light (in the medium), r is the 
distance from the target, and A is the wavelength in the medium (eq.fT2l). See ij24lfor a detailed discussion 
of the Mueller scattering matrix. 

The time-averaged radiative force and torque (see fJT5]l on a target in a dielectric medium are 



F,ad = QprTTa^gUrad , (14) 

Trad = QrTracff^radT;— , (15) 
ZTT 

where the time-averaged energy density is 

El 

^rad — ^medium j (16) 

8tt 

where Eq cos{ujt + (f) is the electric field of the incident plane wave in the medium. 



4 What's New? 

DDSCAT 7.0 differs from DDSCAT 6.1 in a number of ways, including; 

1. N.B.: The structure of the parameter file ddscat .par has been changed: ddscat .par files 
that were used with DDSCAT 6.1 will need to be modified. See Sj8]and AppendixlAl 

2. DDSCAT 7.0 can calculate scattering by targets consisting of anisotropic materials with arbitrary 
orientation. Two new target options ANIFRMFIL and SPH_ANI_N are provided that make use of 
this capability. 

ANIFRMFIL reads an input file with a list of dipoles for a general anisotropic target. For each 
dipole, the input file gives the dipole location, composition for each of three principal axes of 
the local dielectric tensor, and the orientation of the local "Dielectric Frame" (DF) relative to the 
"Target Frame" (TF). 

Target option SPH_ANI_N creates target consisting of N anisotropic spheres. An input file lists, 
for each sphere, the location, radius, composition corresponding to each of three principal axes of 
dielectric tensor, and the orientation of the DF for this sphere relative to the TF. 

3. The user now has the option of specifying the scattering directions either relative to the Lab Frame 
(as in DDSCAT 6. 1) or relative to the Target Frame. This is selected by providing a new character 
variable, either LFRAME or TFRAME in the parameter file ddscat .par. Note that this choice 
must be specified in ddscat . par. 

4. The user must now specify the value of parameter IWRPOL in the parameter file ddscat . par. 
See fj8]and Appendix|A] 

5. A separate Fortran-90 code DDf ield is provided that can read the stored polarization array that 
is written if IWRPOL = 1 and use this data to calculate E and B at user-specified positions near 
the target. 

6. DDSCAT 7.0 can now solve the problem of scattering and absorption by a target that is periodic 
in one dimension (and finite in the other two dimensions) or periodic in two dimensions (and finite 



OBTAINING AND INSTALLING THE SOURCE CODE 



in the third) and illuminated by a monochromatic plane wave, using periodic boundary condi- 
tions (PBC). The theory is described in Draine & Flatau (2008). A number of built-in target op- 
tions are provided (see including BISLINPBC, CYLNDRPBC, DSKLYRPBC, HEXGONPBC, 
LYRSLBPBC, RCTGL_PBC, SPHRN_PBC, 

7. DDSCAT 7.0 now has has dynamic memory allocation. An additional line is now required in 
the ddscat . par file to inform DDSCAT 7.0 how much memory should initially be allocated to 
safely generate the requested target. Once the target has been generated, DDSCAT 7.0 reevaluates 
the memory requirements, and then allocates only the memory required to carry out the scattering 
calculation for the chosen target. 

8. The user can choose to compile DDSCAT 7.0 in either single- or double-precision. 

9. DDSCAT 7.0 includes built-in support for OpenMP; if this is enabled, certain parts of the calcula- 
tion can be distributed over available cpu cores. This can speed execution on multi-core systems. 

10. DDSCAT 7.0 includes built-in support to use the Intel® Math Kernel Library (MKL) for calcula- 
tion of FFTs. 

5 Obtaining and Installing the Source Code 

The complete source code for DDSCAT 7.0 is provided in a single gzipped tarfile. 
There are 3 ways to obtain the source code: 

• Point your browser to the DDSCAT home page: 

|http : / / www . astro . princet on . edu/ ^draine /DDSCAT . html| 
and right=click on the link to ddscat7 .0.7. tgz . 

• Point your browser to 

[ftp : / / ftp . astro .princeton . edu/ draine/ scat/ddscat/7 . . 7| 
right-click on ddscat7 . . 7 . tgz or ddscat7 . . 7 . tgz . 

• Use anonymous ftp to ftp.astro.princeton.edu. 
In response to "Name:" enter anonymous. 

In response to "Password:" enter your full email address. 

cd draine/scat/ddscat/ver7 . . 7 
binary ( to enable transfer of binary data) 
get ddscat7 . . 7 . tgz 

After downloading ddscat7 .0.7. tgz into the directory where you would like DDSCAT 7.0 to re- 
side (you should have at least 10 Mbytes of disk space available), the source code can be installed as 
follows: 

5.1 Installing the source code on a Linux/Unix system 

If you are on a Linux system, you should be able to type 

tar xvzf ddscat7 . . 7 . tgz 
which wiir'extract" the files from the gzipped tarfile. If your version of "tar" doesn't support the "z" 
option (e.g., you are running Solaris) then try 

zcat ddscat7 . . 7 . tgz | tar xvf - 
If neither of the above work on your system, try the two-stage procedure 

gunzip ddscat7 . . 7 . tgz 

tar xvf ddscat7 . . 7 . tar 
A disadvantage of this two-stage procedure is that it uses more disk space, since after the second step 
you will have the uncompressed tarfile ddscat7 .0.7. tar - about 3.8 Mbytes - in addition to all the 
files you have extracted from the tarfile - another 4.6 Mbytes. 

Any of the above approaches should create subdirectories src and doc. The source code will be in 
subdirectory src, and documentation in subdirectory doc. 
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5.2 Installing the source code on a MS Windows system 

If you are running Windows on a PC, you will need the "winzip" program]! winzip should be able 
to "unzip" the gzipped tarfile ddscat 7 . . 7 . tgz and "extract" the various files from it automatically. 

One can create a Linux-like environment under Windows using Cygwin0. Cygwin provides an ex- 
tensive collection of tools which support a Linux-like working environment. The Cygwin DLL currently 
works with recent, comercially-released x86 32 bit and 64 bit versions of windows, and allows users to 
compile and run codes with g77 and g95. 

6 Compiling and Linking 

In the discussion below, it is assumed that the DDSCAT 7.0 source code has been installed in a directory 
DDA/src. The instructions below assume that you are on a Linux system, or in a Cygwin environment 
emulating Linux. If you are not, please see %.^\ 

6.1 The default Makefile 

It is assumed that the system has the following already installed: 

• a Fortran-90 compiler (e.g., gfortran, g95, Intel® ifort,or NAG® f 95) . 

• cpp - the "C preprocessor". 

There are a number of different ways to create an executable, depending on what options the user wants: 

• what compiler and compiler flags? 

• single- or double-precision? 

• enable OpenMP? 

• enable MKL? 

• enable MPl? 

Each of the above choices needs requires setting of appropriate "flags" in the Makefile. 
The default Makefile has the following "vanilla" settings: 

• gfortran -02 

• single-precision arithmetic 

• OpenMP not used 

• MKL not used 

• MPl not used. 

To compile the code with the default settings, simply position yourself in the directory DDA/src, and 
type 

make ddscat 

If you have gfortran and cpp installed on your system, the above should work. You will get some 
warnings from the compiler, but the code will compile and create an executable ddscat. 

If you wish to use a different compiler (or compiler options) you will need to edit the file Make file 
to change the choice of compiler (variable FC), compilation options (variable F FLAGS), and possibly 
and loader options (variable LDFLAGS). The file Makefile as provided includes examples for com- 
pilers other than gfortran; you may be able to simply "comment out" the section of Makefile that 
was designed for gfortran, and "uncomment" a section designed for another compiler (e.g., Intel® 
if ort). 

winzip can be downloaded from |http : / /www . winzip . com| 
|http : / / www . cygwin . com| 
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6.2 Optimization 

The performance of DDSCAT 7.0 will benefit from optimization during compilation and the user should 
enable the appropriate compiler flags. 

6.3 Single vs. Double Precision 

DDSCAT 7.0 is written to allow the user to easily generate either single- or double-precision versions of 
the executable. For most purposes, the single-precision precision version of DDSCAT 7.0 should work 
fine, but for problems where the single-precision version of DDSCAT 7.0 seems to be having trouble 
converging to a solution, the user may wish to experiment with using the double-precision version - this 
can be beneficial in the event that round-off error is compromising the performance of the conjugate- 
gradient solver. Of course, the double precision version will demand about twice as much memory as 
the single-precision version, and will take somewhat longer per iteration. 
The only change required is in the Makefile: for single-precision, set 

PRECISION = sp 
or for double-precision, set 

PRECISION = dp 

After changing the PRECISION variable in the Makefile (either sp — + dp, or dp sp), it is necessary 
to recompile the entire code. Simply type 

make clean 

make ddscat 
to create ddscat with the appropriate precision. 

6.4 OpenMP 

OpenMP is a standard for support of shared-memory parallel programming, and can provide a perfor- 
mance advantage when using DDSCAT 7.0 on platforms with multiple cpus or multiple cores. OpenMP 
is supported by many common compilers, such as gf ortran and Intel® if ort. 

If you are using a multi-cpu (or multi-core) system with OpenMP (www . openmp . org) installed, 
you can compile DDSCAT 7.0 with OpenMP directives to parallelize some of the calculations. To do 
so, simply change 

DOMP = 

OPENMP = 

to 

DOMP = -Dopenmp 
OPENMP = -openmp 
Note: OPENMP is compiler-dependent: gf ortran, for instance, requires 

OPENMP = -fopenmp 



6.5 MKL: the Intel® Math Kernel Library 

Intel® offers the Math Kernel Library (MKL) with the ifort compiler This library includes DFTI for 
computing FFTs. At least on some systems, DFTI offers better performance than the GPFA package. 
To use the MKL library routine DFTI: 

• You must obtain the routine mkl_dfti . f 90 and place a copy in the directory where you are 
compiling DDSCAT 7.0. mkl_df t i . f 9 is Intel® proprietary software, so we cannot distribute 
it with the DDSCAT 7.0 source code, but it should be available on any system with a license for 
the Intel® MKL library. If you cannot find it, ask your system administrator 

• Edit the Makefile: define variables CXFFTMKL . f and CXFFTMKL . o to: 

CXFFTMFKL.f = $(MKL_f) 
CXFFTMKL. f = $(MKL_o) 
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• Successful linking will require that the appropriate MKL libraries be available, and that the string 
LFLAGS in the Makefile be defined so as to include these libraries. 

-Imk:l_em64t -lmkl_intel_thread -lmkl_core -Iguide -Ipthread -lmkl_intel_lp64 
appears to work on at least one installation. The Makefile contains examples. 

• type 

make clean 
make ddscat 

• The parameter file ddscat . par should have FFTMKL as the value of CMETHD. 

6.6 MPI: Message Passing Interface 

DDSCAT 7.0 includes support for parallelization under MPI. MPI (Message Passing Interface) is a 
standard for communication between processes. More than one implementation of MPI exists (e.g., 
mpich and openmpi). MPI support within DDSCAT 7.0 is compliant with the MPl-1.2 and MPl- 
2 standard^ and should be usable under any implementation of MPI that is compatible with those 
standards. 

Many scattering calculations will require multiple orientations of the target relative to the incident 
radiation. For DDSCAT 7.0, such calculations are "embarassingly parallel", because they are carried out 
essentially independently. DDSCAT 7.0 uses MPI so that scattering calculations at a single wavelength 
but for multiple orientations can be carried out in parallel, with the information for different orientations 
gathered together for averaging etc. by the master process. 

If you intend to use DDSCAT 7.0 for only a single orientation, MPI offers no advantage for DDSCAT 
7.0 so you should compile with MPI disabled. However, if you intend to carry out calculations for 
multiple orientations, and would like to do so in parallel over more than one cpu, and you have MPI 
installed on your platform, then you will want to compile DDSCAT 7.0 with MPI enabled. 

To compile with MPI disabled: in the Makefile, set 
DMPI = 

MIP . f = mpi_fake.f90 
MPI.o = mpi_fake.o 
To compile with MPI enabled: in the Makefile, set 

DMPI = -Dmpi 
MIP . f = $ (MPI_f ) 
MPI . o = $ (MPI.o) 

and edit LFLAGS as needed to make sure that you link to the appropriate MPI library (if in doubt, consult 
your systems administrator). The Makefile in the distribution includes some examples, but library 
names and locations are often system-dependent. Please do not direct questions regarding LFLAGS to 
the authors - ask your sys-admin or other experts familiar with your installation. 

Note that the MPl-capable executable can also be used for ordinary serial calculations using a single 
cpu. 

6.7 Device Numbers IDVOUT and IDVERR 

So far as we know, there are only one operating-system-dependent aspect of DDSCAT 7.0: the device 
number to use for "standard output". 

The variables IDVOUT and IDVERR specify device numbers for "running output" and "error mes- 
sages", respectively. Normally these would both be set to the device number for "standard output" (e.g., 
writing to the screen if running interactively). Variables IDVERR are set by DATA statements in the 
"main" program DDSCAT . f 90 and in the output routine WRIMSG (file wrimsg . f 90). The executable 
statement IDVOUT=0 initiahzes IDVOUT to 0. In the as-distributed version of DDSCAT . f 90, the state- 
ment 

OPEN (UNIT=IDVOUT, FILE=CFLLOG) 



" ]http : / /www . mpi-f orum. org/| 
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causes the output to UNIT=IDVOUT to be buffered and written to the file ddscat . log_nnn, where 
nnn=0 00 for the first cpu, 001 for the second cpu, etc. If it is desirable to have this output unbuffered 
for debugging purposes, (so that the output will contain up-to-the-moment information) simply comment 
out this OPEN statement. 

6.8 Installation of DDSCAT 7.0 on non-Linux systems 

DDSCAT 7.0 is written in standard Fortran-90 with a single extension: it uses the Fortran-95 standard 
library call CPU_TIME to obtain timing information. 

It is possible to run DDSCAT on non-Linux systems. If the Unix/Linux "make" utility is not avail- 
able, here in brief is what needs to be accomplished: 

All of the necessary Fortran code to compile and link DDSCAT 7.0 was included in and should have 
been extracted from the gzipped tarfile ddscat? . . tgz. A subset of the files needs to be compiled 
and linked to produce the ddscat executable. 

The main program is in DDSCAT . f 90 ; it calls a number of subroutines. 

Some of the subroutines exist in more than one version. 

• If support for the Intel® Math Kernel Library is not desired, use the file cxf f t3_mkl_f ake . f 90. 
The resulting code will use the Intel® MKL. 

If MKL support is desired (so that option FFTMKL can be used), then use cxf ft3_mkl . f 90 
and mkl_df ti . f 90. The latter is Intel® proprietary, and should already be on you system if you 
have an Intel® MKL license. 

• If MPI support is not desired, use the file mpi_f ake . f 90 instead of mpi_subs . f 90. The 
resulting code will not support MPI. 

If MPI support is required, you will need to install MPI on your system - see ^ 16.61 After doing so, 
you use mpi_subs . f 90 instead of mpi_f ake . f 90. You will also need to compile 

mpi_bchast_char . f 90 
mpi jDcast_cplx . f 90 
mpi jDcast_int . f 90 
mpi_bcast_int2 . f 90 
mpi_bcast_real . f 90 

Prior to compilation, four of the routines 

ddprecision . f 90 
DDSCAT . f 90 
cgcommon . f 90 
eself . f 90 
scat . f 90 

require pre-processing by either cpp or f pp (or equivalent), with the following flags: 

• -Dsp or -Ddp {to enable either single- or double-precision) 

• -Dmp± (to enable MPI) 

• -Dopenmp (to enable OpenMP) 

For example: to compile in single-precision, without MPI, and without OpenMP: 



cpp 


-p 


-traditional- 


-cpp 


-Dsp 


ddprecision . f 90 ddprecision_cpp . f 90 


cpp 


-p 


-traditional- 


-cpp 


-Dsp 


DDSCAT. f 90 DDSCAT.cpp. f 90 


cpp 


-p 


-traditional- 


-cpp 


-Dsp 


cgcommon. f 90 cgcommon_cpp . f 90 


cpp 


-p 


-traditional- 


-cpp 


-Dsp 


eself. f90 eself _cpp . f 90 


cpp 


-p 


-traditional- 


-cpp 


-Dsp 


scat. f 90 scat_cpp.f90 



Before compiling any other routines, you need to first compile 

ddprecision . f 90 
which will create the module 

ddprecision . mod 
Next compile 

ddcommon . f 90 
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which will create the following modules: 

ddcommon_l .mod 

ddcommon_2 .mod 

ddcommon_3 . mod 

ddcommon_4 .mod 

ddcommon_5 . mod 

ddcommon_6 . mod 

ddcommon_7 .mod 

ddcommon_8 .mod 

ddcommon_9 . mod 

ddcommon_10 .mod 
After these modules have been created, proceed to compile 

DDSCAT_cpp . f 90 

cgcommon_cpp . f 90 

eself_cpp . f 90 

scat_cpp . f 90 

alphadiag . f 90 

blas.f90 

Note: do «of compile the original ddprecision . f 90, DDSCAT . f 90, cgcommon .f90,eself.f90, 
or scat . f 90 - these are superceded by DDSCAT_cpp . f 90, cgcommon_cpp . f 90, eself _cpp . f 90, 
and scat_cpp . f 90. 

Try to link the files DDSCAT_cpp . o and the other . o files to create an executable, just as you would 
with any other Fortran code with a number of modules. You should end up with an executable with a 
(system-dependent) name like DDSCAT_cpp . EXE. You may wish to rename this to, e.g., DDSCAT . EXE. 

If you do not have either cpp or f pp available, you can use an editor to search through the source 
code for 

#if def option 

#endif 

and enable the appropriate blocks of code (and comment out the undesired blocks) and then compile and 
link. 

NOTE: If you have trouble with the above general directions, you are strongly encouraged to install 
Cygwin (see ^15. 21 ). which will allow you to instead use the Makefile as described above for installation 
on Linux systems. 

7 Moving the Executable 

Now reposition yourself into the directory DDA (e.g., type cd . . ), and copy the executable from 
src/ddscat to the DDA directory by typing 

cp src/ddscat ddscat 
This should copy the file DDA/src/ddscat to DDA/ ddscat (you could, alternatively, create a sym- 
bolic link). Similarly, copy the sample parameter file ddscat . par and the file diel . tab to the DDA 
directory by typing 

cp doc/ddscat . par ddscat. par 

cp doc/diel.tab diel. tab 

8 The Parameter File ddscat. par 

The directory DDA should contain a sample file ddscat .par (see Appendix lAli which provides pa- 
rameters to the program ddscat. Appendix lAl provides a line-by-line description of the sample file 
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ddscat . par|l Here we discuss the general structure of ddscat . par. 

8.1 Preliminaries 

ddscat . par starts by setting the values of five strings: 

• CMDTRQ specifying whether or not radiative torques are to be calculated (e.g., NOTORQ) 

• CMDSOL specifying the solution method (e.g., PBCGS2) (see section ^TTT i. 

• CMDFFT specifying the FFT method (e.g., GPFAFT) (see section CT . 

• CALPHA specifying the polarizability prescription (e.g., GKDLDR) (see 23- 

• CBINFLAG specifying whether to write out binary files (e.g., NOTBIN) 

8.2 Initial Memory Allocation 

Three integers are given that need to be equal to or larger than the anticipated target size (in lattice units) 
in the xtf, Ytf, and ztf directions. This is required for initial memory allocation. A value like 
100 100 100 

would require initial memory allocation of only ~ 100 MBytes. Note that after the target geometry has 
been allocated, DDSCAT 7.0 will proceed to reallocate only as much memory as is actually required for 
the calculation. 

8.3 Target Geometry and Composition 

• CSHAPE specifies the target geometry (e.g., RCTGLPRSM) 

As provided (see Appendi?jA]l. the file ddscat . par is set up to calculate scattering by a 32x24x 16 
rectangular array of 12288 dipoles. 

The user must specify NCOMP = how many different dielectric functions will be used. This is then 
followed by NCOMP lines, with each line giving the name (in quotes) of each dielectric function file. In 
the example, the dielectric function of the target material is provided in the file diel . tab, which is a 
sample file in which the refractive index is set to m = 1.33 + O.Oli at all wavelengths. 

8.4 Error Tolerance 

TOL = the error tolerance. Conjugate gradient iteration will proceed until the linear equations are solved 
to a fractional error TOL. The sample calculation has TOL=10~^. 

8.5 Interaction Cutoff Parameter for PBC Calculations 

GAMMA = parameter limiting certain summations that are required for periodic targets (see Draine & 
Flatau 2008). The value of GAMMA has no effect on computations for finite targets - it can be set to 
any value, including 0). 

For targets that are periodic in 1 or 2 dimensions, GAMMA needs to be small for high accuracy, but 
the required cpu time increases as GAMMA becomes smaller. GAMMA = 10^^ is reasonable for initial 
calculations, but you may want to experiment to see if the results you are interested in are sensitive to 
the value of GAMMA. Draine & Flatau (2008) show examples of how computed results for scattering can 
depend on the value of 7. 



* Note that the structure of dds cat .par depends on the version of DDSCAT being used. Make sure you update old parameter 
files before using them with DDSCAT 7.0 ! 
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8.6 Angular Resolution for Computation of (cos 6), etc. 

The parameter ETASCA determines the selection of scattering angles used for computation of certain 
angular averages, such as {cos 9) and {cos^ 9), and the radiation pressure force (see ^031 ) and radia- 
tive torque (if CMDTRQ=DOTORQ). Small values of ETASCA result in increased accuracy but also cost 
additional computation time. ETASCA=0.5 generally gives accurate results. 

If accurate computation of (cos 6) or the radiation pressure force is not required, the user can set 
ETASCA to some large number, e.g. 10, to minimize unnecessary computation. 

8.7 Wavelength 

The sample ddscat . par file specifies that the calculations be done for a single wavelength (6.2832). 
[If one wishes to specify the "size parameter" x = 27racff / A, it is convenient to set A = 27r so that 

X = Oeff.] 

8.8 Target Size Oeff 

Note that in DDSCAT the "effective radius" Ocft is the radius of a sphere of equal volume - i.e., a sphere 
of volume Nd^ , where d is the lattice spacing and N is the number of occupied (i.e., non-vacuum) 
lattice sites in the target. Thus the effective radius Ocff = (3A^/47r)^/'^d . The sample ddscat . par file 
specifies an effective radius — 2/im. If the rectangular solid is a x 6 x c, with a : & : c :: 32 : 24 : 16, 

then abc = 3c^ = (47r/3)a3ff, so that c = {ATT/9)^/^acs = 2.235/im, b = (3/2)c = 3.353^im, 
a = 2c = 4.470^m. 

The "size parameter" 27raoff/A = 2 for the sample ddscat . par. 

8.9 Incident Polarization 

The incident radiation is always assumed to propagate along the x axis in the "Lab Frame". The sample 
ddscat .par file specifies incident polarization state Gqi to be along the y axis (and consequently 
polarization state eo2 will automatically be taken to be along the z axis). I0RTH=2 in ddscat .par 
calls for calculations to be carried out for both incident polarization states (eoi and eo2 - see f|22]i. 

8.10 Output Files to Write 

The sample ddscat .par file sets 1 = IWRKSC, so that ".sea" files will be written out with scattering 
information for each orientation. 

The sample ddscat .par file sets 1 = IWRP0L, so that a ".pol" file will be writtten out for each 
incident polarization and each orientation. Note that each such file has a size that is proportional to the 
number of dipoles, and can be quite large. For this modest test problem, with N = 12288 dipoles, the 
.pol files are each 958 kBytes. These files are only useful if post-processing of the solution is desired, 
such as evaluation of the electromagnetic field near the target using the DDf ield program as described 
ingZl 

8.11 Target Orientation 

The target is assumed to have two vectors ai and 3.2 embedded in it; a2 is perpendicular to ai . In the 
case of the 32 x 24 x 1 6 rectangular array of the sample calculation, the vector ai is along the "long" axis 
of the target, and the vector 3.2 is along the "intermediate" axis. The target orientation in the Lab Frame 
is set by three angles: /3, 8, and <&, defined and discussed below in fJTll Briefly, the polar angles 8 and 
$ specify the direction of ai in the Lab Frame. The target is assumed to be rotated around ai by an angle 
f3. The sample ddscat .par file specifies (3 = and $ = (see lines in ddscat .par specifying 
variables BETA and PHI), and calls for three values of the angle 8 (see line in ddscat . par specifying 
variable THETA). DDSCAT chooses 8 values uniformly spaced in cos 8; thus, asking for three values 
of 8 between and 90° yields 8 = 0, 60°, and 90°. 
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8.12 Starting Values of IWAV, IRAD, lORl 

Normally we begin the calculation with IWAVE=0, IRAD=0, and IORI=0. However, under some cir- 
cumstances, a prior calculation may have completed some of the cases. If so, the user can specify starting 
values of IWAV, IRAD, lORI; the computations will begin with this case, and continue. 

8.13 Which Mueller Matrix Elements? 

The sample parameter file specifies that DDSCAT should calculate 6 distinct Mueller scattering matrix 
elements Sij, with 11, 12, 21, 22, 31, 41 being the chosen values of ij. 

8.14 What Scattering Directions? 

8.14.1 Isolated Finite Targets 

For finite targets, the user may specify the scattering directionsin either the Lab Frame (' LFRAME' ) or 
the Target Frame (' TFRAME' )■ 

For finite targets, such as specified in this sample ddscat .par, the Sij are to be calculated for 
scattering directions specified by angles {6, (f>). 

The sample ddscat . par specifies that NPLANES=2 scattering planes are to be specified: the first 
has <j> = and the second has cj) ~ 90°; for each scattering plane 6 values run from to 180° in 
increments of 10°. 

8.14.2 1-D Periodic Targets 

For periodic targets, the scattering directions must be specified in the Target Frame: ' TFRAME' = 
CMDFRM . 

Scattering from 1-d periodic targets is discussed in detail by Draine & Flatau (2008). For periodic 
targets, the user does not specify scattering planes. For 1 -dimensional targets, the user specifies scattering 
cones, corresponding to different scattering orders M. For each scattering cone M, the user specifies 
Cmin, Cmax, (i^ degrees); the azimuthal angle ( will run from (jnin to C,rnax, in increments of A^. For 
example: 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 
1 = NPLANES = number of scattering cones 

0. 0. 180. 0.05 = OrderM zetamin zetamax dzeta for scattering cone 1 

8.14.3 2-D Periodic Targets 

For targets that are periodic in 2 dimensions, the scattering directions must be specified in the Target 
Frame: ' TFRAME' = CMDFRM. 

Scattering from targets that are periodic in 2 dimensions is discussed in detail by Draine & Flatau 
(2008). For 2-D periodic targets, the user specifies the diffraction orders {M, N) for transmitted radiation 
(the code will automatically also calculate reflection coefficients for the same scattering orders.) For 
example: 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = NPLANES = number of scattering cones 

0. 0. = OrderM OrderN for scattered radiation 

9 Running DDSCAT 7.0 Using the Sample dd scat. par File 
9.1 Single-Process Execution 

To execute the program on a UNIX system (running either sh or csh), simply type 
./ddscat >& ddscat. out & 
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which will redirect the "standard output" to the file ddscat . out, and run the calculation in the back- 
ground. The sample calculation [32x24x16=12288 dipole target, 3 orientations, two incident polar- 
izations, with scattering (Mueller matrix elements Sij) calculated for 38 distinct scattering directions], 
requires 27.9 cpu seconds to complete on a 2.0 GHz single-core Xeon workstation. 

9.2 Code Execution Under MPI 

Local installations of MP I will vary - you should consult with someone familiar with way MP I is in- 
stalled and used on your system. 

At Princeton University Dept. of Astrophysical Sciences we use PBS (Portable Batch Systemj^l to 
schedule jobs. MPI jobs are submitted using PBS by first creating a shell script such as the following 
example file pbs . submit: 

#!/bin/bash 

#PBS -1 nodes=2:ppn=l 

#PBS -1 mem=1200MB,pmem=300MB 

#PBS -m bea 

#PBS -j oe 

cd $PBS_0_WORKDIR 

/usr/local/bin/mpiexec ddscat 

The Unes beginning with #PBS -1 specify the required resources: 

#PBS -1 nodes=2 : ppn=l specifies that 2 nodes are to be used, with 1 processor per node. 
#PBS -1 mein=12 00MB, pinem=300MB specifies that the total memory required (mem) is 1200MB, 
and the maximum physical memory used by any single process (pmem) is 300MB. The actual definition 
of mem is not clear, but in practice it seems that it should be set equal to 2 x (node s) x (ppn) x (pmem) . 
#PBS -m bea specifies that PBS should send email when the job begins (b), and when it ends (e) or 
aborts (a). 

#PBS - j oe specifies that the output from stdout and stderr will be merged, intermixed, as stdout. 
This example assumes that the executable dscat is located in the same directory where the code is to 
execute and write its output. If ddscat is located in another directory, simply give the full pathname, 
to it. The qsub command is used to submit the PBS job: 

qsub pbs. submit 

As the calculation proceeds, the usual output files will be written to this directory: for each wave- 
length, target size, and target orientation, there will be a file waaarbbbkccc . sea, where aaa=000, 001, 
002, ... specifies the wavelength, bbb=000, 001, 002, ... specifies the target size, and cc=000, 001, 002, ... 
specifies the orientation. For each wavelength and target size there will also be a file waaarbbbori . avg 
with orientationally-averaged quantities. Finally, there will also be tables qtable, and qtable2 with 
orientationally-averaged cross sections for each wavelength and target size. 

In addition, each processor employed will write to its own log file dds cat . log jtnn, where nnn=000, 
001, 002, ... . These files contain information concerning cpu time consumed by different parts of the cal- 
culation, convergence to the specified error tolerance, etc. If you are uncertain about how the calculation 
proceeded, examination of these log files is recommended. 

10 Output Files 
10.1 ASCII files 

If you run DDSCAT using the command 

ddscat >& ddscat. out & 

' ^http : / / www ■ openpbs . org| 
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you will have various types of ASCII files when the computation is complete: 

• a file ddscat . out; 

• a file ddscat . log_000; 

• a file mtable; 

• a file qtable; 

• a file qtable2; 

• files vxxxryyori . avg (one, wOOOrOOOori . avg, for the sample calculation); 

• if ddscat .par specified IWRKSC=1, there will also be files wxxxryykzzz . sea (3 for the sample 
calculation: wOOOrOOOkOOO . sea, wOOOrOOOkOOl . sea, w000r000k002 . sea). 

• if ddscat .par specified IWRP0L=1, there will also be files vixxxryy^kziz . poln forn=l and (if 

I0RTH=2) n=2 

The file ddscat . out will contain minimal information (it may in fact be empty). 

The file ddscat . log_000 will contain any error messages generated as well as a running report 
on the progress of the calculation, including creation of the target dipole array. During the iterative 
calculations, Qext, Qabs, and Qpha are printed after each iteration; you will be able to judge the degree 
to which convergence has been achieved. Unless TIME IT has been disabled, there will also be timing 
information. If the MP I option is used to run the code on multiple cpus, there will be one file of the form 
ddscat . log_nnn for each of the cpus, with nnn=000 , 001,002, . . .. 

The file mtable contains a summary of the dielectric constant used in the calculations. 

The file qtable contains a summary of the orientationally-averaged values of Qoxt, Qabs, Qsca, 
g{l) ^ {cos{8s)), (cos2(6's)), Qbk, and iVsca- Here Qext, Qabs, and Qsca are the extinction, absorption, 
and scattering cross sections divided by Tra^jj . Qbk is the differential cross section for backscattering 
(area per sr) divided by vra^fj. A^sca is the number of scattering directions used for averaging over 
scattering directions (to obtain (cos 9), etc.) (see ^23^ . 

The file qtable2 contains a summary of the orientationally-averaged values of Qpha, Qpoh and 
Qcpoi- Here Qpha is the "phase shift" cross section divided by Tra^fj (see definition in Draine 1988). 
Qpoi is the "polarization efficiency factor", equal to the difference between Qcxt for the two orthogonal 
polarization states. We define a "circular polarization efficiency factor" Qcpoi = QpoiQpha, since an 
optically-thin medium with a small twist in the alignment direction will produce circular polarization in 
initially unpolarized light in proportion to Qcpoi- 

For each wavelength and size, DDSCAT 7.0 produces a file with a name of the form 
vixxxryyori . avg, where index xxx (=000, 001, 002....) designates the wavelength and index yy (=00, 
01, 02...) designates the "radius"; this file contains Q values and scattering information averaged over 
however many target orientations have been specified (see jfTTl The file wOOOrOOOori.avg produced 
by the sample calculation is provided below in AppendixlB] 

In addition, if ddscat .par has specified IWRKSC=1 (as for the sample calculation), DDSCAT 
7.0 will generate files with names of the form wxxxryykzzz . avg, where xxx and yy are as before, and 
index zzz =(000,001,002...) designates the target orientation; these files contain Q values and scattering 
information for each of the target orientations. The structure of each of these files is very similar to 
that of the vixxxryyori . avg files. Because these files may not be of particular interest, and take up 
disk space, you may choose to set IWRKSC=0 in future work. However, it is suggested that you run the 
sample calculation with IWRKSC=1. 

The sample ddscat . par file specifies IWRKSC=1 and calls for use of 1 wavelength, 1 target size, 
and averaging over 3 target orientations. Running DDSCAT 7.0 with the sample ddscat . par file will 
therefore generate files wOOOrOOOkOOO . sea, wOOOrOOOkOOl . sea, and w000r000k002 . sea . 
To understand the information contained in one of these files, please consult AppendixICl which contains 
an example of the file wOOOrOOOkOOO. sea produced in the sample calculation. 

10.2 Binary Option 

It is possible to output an "unformatted" or "binary" file (dd.bin) with fairly complete information, 
including header and data sections. This is accomplished by specifying either ALLBIN or ORIBIN in 
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ddscat . par . 

Subroutine writebin . f 90 provides an example of how this can be done. The "header" section 
contains dimensioning and other variables which do not change with wavelength, particle geometry, 
and target orientation. The header section contains data defining the particle shape, wavelengths, parti- 
cle sizes, and target orientations. If ALLBIN has been specified, the "data" section contains, for each 
orientation, Mueller matrix results for each scattering direction. The data output is limited to actual di- 
mensions of arrays; e.g. nscat, 4, 4 elements of Muller matrix are written rather than mxs cat , 4, 4. 
This is an important consideration when writing postprocessing codes. 



11 Choice of Iterative Algorithm 

As discussed elsewhere (e.g., Draine 1988), the problem of electromagnetic scattering of an incident 
wave Einc by an array of point dipoles can be cast in the form 

AP = E (17) 

where E is a 3A^-dimensional (complex) vector of the incident electric field Ejnc at the N lattice sites, 
P is a 3A^-dimensional (complex) vector of the (unknown) dipole polarizations, and A is a 3N x 3A'' 
complex matrix. 

Because 3N is a large number, direct methods for solving this system of equations for the unknown 
vector P are impractical, but iterative methods are useful: we begin with a guess (typically, P = 0) for 
the unknown polarization vector, and then iteratively improve the estimate for P until equation (fTTI l is 
solved to some error criterion. The error tolerance may be specified as 

|AtAP-AtE| 

\A^\ <^ ' ^^^^ 

where A^ is the Hermitian conjugate of A [(^^^y = (Aji)*], and h is the error tolerance. We typically 
use h = 10^^ in order to satisfy eq.dTTb to high accuracy. The error tolerance h can be specified by the 
user through the parameter TOL in the parameter file ddscat . par (see AppendixlAli. 

A major change in going from DDSCAT.4b to 5a (and subsequent versions) was the implementation 
of several different algorithms for iterative solution of the system of complex linear equations. DDSCAT 
7.0 is now structured to permit solution algorithms to be treated in a fairly "modular" fashion, facilitating 
the testing of different algorithms. A number of algorithms were compared by Flatau two of 

them (PBCGST and PETRKP) performed well and are made available to the user in DDSCAT. DDSCAT 
7.0 now also offers a third option, PBCGS2. The choice of algorithm is made by specifying one of the 
options: 

• PBCGS2 - BiConjugate Gradient with Stabilization as implemented in the routine ZBCG2 by 
M.A. Botchev, University of Twente. This is based on the PhD thesis of D.R. Fokkema, and on 
work by Sleijpen & vanderVorst (1995,1996). 

• PBCGST - Preconditioned BiConjugate Gradient with STabihtization method from the Parallel 
Iterative Methods (PIM) package created by R. Dias da Cunha and T. Hopkins. 

• PETRKP - the complex conjugate gradient algorithm of Petravic & Kuo-Petravic (1979), as coded 
in the Complex Conjugate Gradient package (CCGPACK) created by P.J. Flatau. This is the algo- 
rithm discussed by Draine (1988) and used in previous versions of DDSCAT. 

All three methods work fairly well. Our experience suggests that PBCGS2 is generally fastest and best- 
behaved, and we recommend that the user try it first. However, in case it runs into numerical difficulties, 
PBCGST and PETRKP are available as alternatives0 

^ A postscript copy of this report - file eg . ps - is distributed with the DDSCAT 7.0 documentation. 

' The Parallel Iterative Methods (PIM) by Rudnei Dias da Cunha (rddSukc . ac . uk) and Tim Hopkins (trhSukc . ac . uk) 
is a collection of Fortran 77 routines designed to solve systems of linear equations on parallel and scalar computers using a variety 
of iterative methods (available at 
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12 Choice of FFT Algorithm 

DDSCAT 7.0 offers two FFT options: (1) the GPFA FFT algorithm developed by Dr Clive Temperton 
(1992),0and (2) the Intel® MKL routine DFTI. 

The GPFA routine is portable and quite fast: Figure 5 compares the speed of three FFT implemen- 
tations: Brenner's, GPFA, and FFTW (http : //www . f ftw . org). We see that while for some cases 
FFTW 2 . 1 . 5 is faster than the GPFA algorithm, the difference is only marginal. The FFTW code and 
GPFA code are quite comparable in performance - for some cases the GPFA code is faster, for other 
cases the FFTW code is faster. For target dimensions which are factorizable as 2*3^5*^ (for integer i, 
j, k), the GPFA and FFTW codes have the same memory requirements. For targets with extents N^, 
Ny, Nz which are not factorizable as 2*3^5*^, the GPFA code needs to "extend" the computational vol- 
ume to have values of N^, Ny, and which are factorizable by 2, 3, and 5. For these cases, GPFA 
requires somewhat more memory than FFTW. However, the fractional difference in required memory is 
not large, since integers factorizable as 2*3^ 5*^ occur fairly frequently03 [Note: This "extension" of the 
target volume occurs automatically and is transparent to the user.] 

DDSCAT 7.0 offers a new FFT option: the Intel® Math Kernel Library DFTI. This is tuned for 
optimum performance, and appears to offer real performance advantages on modern multi-core cpus. 
With this now available, the FFTW option, which had been included in DDSCAT 6.1, has been removed 
from DDSCAT 7.0. 

The choice of FFT implementation is obtained by specifying one of: 

|http : / /www ■ mat ■ uf rgs . br /pim-e . html| >. PIM offers a number of iterative methods, including 

• Conjugate-Gradients, CG (Hestenes 1952), 

• Bi-Conjugate-Gradients, BIGG (Fletcher 1976), 

• Conjugate-Gradients squared, CGS (Sonneveld 1989), 

• the stabilised version of Bi-Conjugate-Gradients, BICGSTAB (van der Vorst 1991), 

• the restarted version of BICGSTAB, RBICGSTAB (Sleijpen & Fokkema 1992) 

• the restarted generalized minimal residual, RGMRES (Saad 1986), 

• the restarted generalized conjugate residual, RGCR (Eisenstat 1983), 

• the normal equation solvers, CGNR (Hestenes 1952 and CGNE (Craig 1955), 

• the quasi-minimal residual, QMR (highly parallel version due to Bucker & Sauren 1996), 

• transpose-free quasi-minimal residual, TFQMR (Freund 1992), 

• the Chebyshev acceleration, CHEBYSHEV (Young 1981). 

The source code for these methods is distributed with DDSCAT but only PBCGST and PETRKP can be called directly via 
ddscat.par. It is possible to add other options by changing the code in getfml.f90 . Flatau (1997) has compared 
the convergence rates of a number of different methods. A helpful introduction to conjugate gradient methods is provided 
by the report "Conjugate Gradient Method Without Agonizing Pain" by Jonathan R. Shewchuk, available as a postscript file: 
|ftp : //REPORTS .ADM. CS ■ CMU . EDU/usr /anon/ 1 9 94 /CMU-CS-94-1 2 5 ■ ps| 

The GPFA code contains a parameter LVR which is set in data statements in the routines gpf a2f , gpf a3f , and gpf a5f . 
LVR is supposed to be optimized to correspond to the "length of a vector register" on vector machines. As delivered, this parameter 
is set to 64, which is supposed to be appropriate for Grays other than the C90. For the C90, 128 is supposed to be preferable (and 
perhaps "preferable" should be read as "necessary" - there is some basis for fearing that results computed on a C90 with LVR other 
than 128 run the risk of being incorrect!) The value of LVR is not critical for scalar machines, as long as it is fairly large. We found 
little difference between LVR=64 and 128 on a Sparc 10/51, on an Ultrasparc 170, and on an Intel® Xeon cpu. You may wish to 
experiment with different LVR values on your computer architecture. To change LVR, you need to edit gpf a . f 9 and change the 
three data statements where LVR is set. 

"2, 3, 4, 5, 6, 8, 9, 10, 12, 15, 16, 18, 20, 24, 25, 27, 30, 32, 36, 40, 45, 48, 50, 54, 60, 64, 72, 75, 80, 81, 90, 96, 100, 108, 120, 
125, 128, 135, 144, 150, 160, 162, 180, 192, 200 216, 225, 240, 243, 250, 256, 270, 288, 300, 320, 324, 360, 375, 384, 400, 405, 
432, 450, 480, 486, 500, 512, 540, 576, 600, 625, 640, 648, 675, 720, 729, 750, 768, 800, 810, 864, 900, 960, 972, 1000, 1024, 
1080, 1125, 1152, 1200, 1215, 1250, 1280, 1296, 1350, 1440, 1458, 1500, 1536, 1600, 1620, 1728, 1800, 1875, 1920, 1944, 2000, 
2025, 2048, 2160, 2187, 2250, 2304, 2400, 2430, 2500, 2560, 2592, 2700, 2880, 2916, 3000, 3072, 3125, 3200, 3240, 3375, 3456, 
3600, 3645, 3750, 3840, 3888, 4000, 4050, 4096 are the integers < 4096 which are of the form 2'3-' 5'=. 
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Figure 5: Comparison of cpu time required by 3 different FFT implementations. It is seen that the GPFA and FPTW implemen- 
tations have comparable speeds, much faster than Brenner's FFT implementation. 



• FFTMKL to use the Intel® MKL routine DFT I (see this is recommended, but requires that 
the Intel® Math Kernel Library be installed on your system. 

• GPFAFT to use the GPFA algorithm (Temperton 1992) - this is is not as fast as FFTMKL, but 
is written in plain Fortran-90. 



13 Dipole Polarizabilities 

Option GKDLDR specifies that the polarizability be prescribed by the "Lattice Dispersion Relation", 
with the polarizability found by Gutkowicz-Krusin & Draine (2004), who corrected a subtle error in the 
analysis of Draine & Goodman (1993). For |TO|fcd ^ 1, the GKDLDR polarizability differs slightly from 
the LATTDR polarizability, but the differences in calculated scattering cross sections are relatively small, 
as can be seen from Figure|6] We recommend option GKDLDR. 

Users wishing to compare can invoke option LATTDR to specify that the "Lattice Dispersion Rela- 
tion" of Draine & Goodman (1993) be employed to determine the dipole polarizabilities. This polariz- 
ability also works well. 



14 Dielectric Functions 

In order to assign the appropriate dipole poIarizabiUties, DDSCAT 7.0 must be given the dielectric 
constant of the material (or materials) of which the target of interest is composed relative to the dielectric 
constant of the ambient medium. This information is supplied to DDSCAT 7.0 through a table (or tables), 
read by subroutine DIELEC in file dielec.f90, and providing either the complex refractive index 
m ~ n + ikor complex dielectric function e = €i+ it2 as a function of wavelength A. Since m = e^/^, 
or e = in?, the user must supply either m or e. 

DDSCAT 7.0 can calculate scattering and absorption by targets with anisotropic dielectric functions, 
with arbitrary orientation of the optical axes relative to the target shape. See i j26] 
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Figure 6: Scattering and absorption for a m = 1.7 + O.li sphere, calculated using two prescriptions for the polarizability: 
LATTDR is the lattice dispersion relation result of Draine & Goodman (1993). GKDLDR is the lattice dispersion relation result 
of Gutkowicz-Ki'usin & Draine (2004). Results are shown as a function of scattering parameter x — 27ra/A; the upper scale gives 
values of |m|fcd. We see that the cross sections calculated with these two prescriptions are quite similar for \m\kd ^ 0.5. For other 
examples see Gutkowicz-Krusin & Draine (2004). 



As discussed in DDSCAT calculates scattering for targets embedded in dielectric media us- 
ing the relative dielectric constant e = etargot(ij-')/emcdium(w) or the relative refractive index m = 
mtarget(w)/minedium(w) and the Wavelength in the medium A = (27rc/ct')/m„iedium, where c is the 
speed of light /« vacuo. 

It is therefore important that the table containing the dielectric function information give the relative 
dielectric function or refractive index as a function of the wavelength in the ambient medium. 

The table formatting is intended to be quite flexible. The first line of the table consists of text, up 
to 80 characters of which will be read and included in the output to identify the choice of dielectric 
function. (For the sample problem, it consists of simply the statement m = 1.33 + O.Oli.) The 
second line consists of 5 integers; either the second and third or the fourth and fifth should be zero. 

• The first integer specifies which column the wavelength is stored in. 

• The second integer specifies which column Re(m) is stored in. 

• The third integer specifies which column Im(TO) is stored in. 

• The fourth integer specifies which column Re(e) is stored in. 

• The fifth integer specifies which column Im(e) is stored in. 

If the second and third integers are zeros, then DIELEC will read Re(e) and Im(e) from the file; if the 
fourth and fifth integers are zeros, then Re(m) and Im(m) will be read from the file. 

The third line of the file is used for column headers, and the data begins in line 4. There must be at 
least 3 lines of data: even if e is required at only one wavelength, please supply two additional "dummy" 
wavelength entries in the table so that the interpolation apparatus will not be confused. 

15 Calculation of (cos 6), Radiative Force, and Radiation Torque 

In addition to solving the scattering problem for a dipole array, DDSCAT can compute the three- 
dimensional force Fiad and torque r,ad exerted on this array by the incident and scattered radiation 
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fields. The radiation torque calculation is carried out, after solving the scattering problem, only if 
DOTORQ has been specified in ddscat .par. For each incident polarization mode, the results are 
given in terms of dimensionless efficiency vectors Qp,. and Qr, defined by 

Qpr ^ I"'' , (19) 

Q,^ , (20) 

where Frad and Fiad are the time-averaged force and torque on the dipole array, k = 27r/A is the 
wavenumber in vacuo, and Uiad = Eq/8tt is the time-averaged energy density for an incident plane 
wave with amplitude Eq cos{ujt + (f). The radiation pressure efficiency vector can be written 

Qpr = Qoxtk — Qscag i (21) 

where k is the direction of propagation of the incident radiation, and the vector g is the mean direction 
of propagation of the scattered radiation: 

where dVl is the element of solid angle in scattering direction n, and dCsc&/d^ is the differential scat- 
tering cross section. The components of Qpi are reported in the Target Frame: Qpi.i ^ Fiad • xtf, 

Qpr, 2 = Frad ' YTF, Qpr, 3 = Frad ' ^TF- 

Equations for the evaluation of the radiative force and torque are derived by Draine & Weingartner 
(1996). It is important to note that evaluation of Qpr and Qr involves averaging over scattering directions 
to evaluate the linear and angular momentum transport by the scattered wave. This evaluation requires 
appropriate choices of the parameter ETASCA - see f|23] 

In addition, DDSCAT calculates (cos 0) [the first component of the vector g in eq. (|22]) 1 and the 
second moment (cos^ 0) . These two moments are useful measures of the anisotropy of the scattering. 
For example, Draine (2003) gives an analytic approximation to the scattering phase function of dust 
mixtures that is parameterized by the two moments (cos 9) and (cos^ 0) . 



16 Memory Requirements 

The memory requirements are determined by the size of the "computational volume" - this is a rectan- 
gular region, of size NX x NY x NZ that is large enough to contain the target. If using the GPFAFT option, 
then NX, NY, NZ are also required to have only 2,3, and 5 as prime factors (see footnote [TTli. 
In single precision, the memory requirement for DDSCAT 7.0 is approximately 

(35. + 0.0010 X NX X NY X NZ) Mbytes for single precision (23) 

(42 + 0.0020 X NX X NY X NZ) Mbytes for double precision (24) 

Thus, in single precision, a 48x48x48 calculation requires ~146 MBytes. 

The memory is allocated dynamically - once the target has been created, DDSCAT 7.0 will deter- 
mine just how much overall memory is needed, and will allocate it. However, the user must provide 
information (via ddscat . par) to allow DDSCAT 7.0 to allocate sufficiently large arrays to carry out 
the initial target creation. Initially, the only arrays that will be allocated are those related to the target 
geometry, so it is OK to be quite generous in this initial allowance, as the memory required for the target 
generation step is small compared to the memory required to carry out the full scattering calculation. 
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Fi gure 7: Target orientation in the Lab Frame, x = xlf is the direction of propagation of the incident radiation, and y — ylf 
is the direction of the real component (at slf ~ 0, t — 0) of the first incident polarization mode. In this coordinate system, the 
orientation of target axis ai is specified by angles B and With target axis ai fixed, the orientation of target axis a2 is then 
determined by angle P specifying rotation of the target around ai. When /3 = 0, a2 lies in the ai,XLF plane. 



17 Target Orientation 

Recall that we define a "Lab Frame" (LF) in which the incident radiation propagates in the +x direction. 
For purposes of discussion we will always let unit vectors xlf, Ylf, zlf = xlf x Ylf be the three 
coordinate axes of the LF. 

In ddscat .par one specifies the first polarization state Gqi (which obviously must lie in the y, z 
plane in the LF); DDSCAT automatically constructs a second polarization state eo2 = xlf x ^oi orthog- 
onal to Gqi. Users will often find it convenient to let polarization vectors epi = y, eo2 = z (although 
this is not mandatory - see i j22] i. 

Recall that definition of a target involves specifying two unit vectors, ai and 3.2, which are imagined 
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to be "frozen" into the target. We require a.2 to be orthogonal to ai . Therefore we may define a "Target 
Frame" (TF) defined by the three unit vectors ai, a2, and as = ai x 3.2 . 

For example, when DDSCAT creates a 32x24x16 rectangular solid, it fixes ai to be along the 
longest dimension of the solid, and a2 to be along the next-longest dimension. 

Important Note: for periodic targets, DDSCAT 7.0 requires that the periodic target have ai = xtf 
and a2 = yxF ■ 

Orientation of the target relative to the incident radiation can in principle be determined two ways: 

1 . specifying the direction of ai and a2 in the LF, or 

2. specifying the directions of xlf (incidence direction) and ylf in the TF. 

DDSCAT uses method 1.: the angles O, $, and /3 are specified in the file ddscat .par. The target is 
oriented such that the polar angles 8 and $ specify the direction of ai relative to the incident direction 
Xlf; where the xlf.Ylf plane has $ = 0. Once the direction of ai is specified, the angle /3 then 
specifies how the target is to rotated around the axis ai to fully specify its orientation. A more extended 
and precise explanation follows: 



17.1 Orientation of the Target in the Lab Frame 

DDSCAT uses three angles, 8, <&, and f3, to specify the directions of unit vectors ai and a2 in the LF 
(see Fig.|7|i. 

8 is the angle between ai and xlf- 

When <1) = 0, ai will lie in the xlf, Ylf plane. When $ is nonzero, it will refer to the rotation of ai 
around xlf^ e.g., $ = 90° puts ai in the xlf, zlf plane. 

When /? = 0, 3.2 will lie in the xlf , ai plane, in such a way that when 8 = and $ = 0, a2 is in 
the yLF direction: e.g, 8 ~ 90°, $ = 0, /3 = has ai jlf and 3.2 = — xlf- Nonzero [3 introduces 
an additional rotation of a2 around ai: e.g., 8 = 90°, <i> = 0, /? = 90° has ai = ylf and 3.2 = zlf- 

Mathematically: 

ai = Xlf cos 8 + ylf sin 8 cos $ + zlf sin 8 sin $ (25) 
62 = —Xlf sin 8 cos (3 + ylf [cos 8 cos (3 cos $ — sin /3 sin $] 

+ZLF [cos 8 cos (3 sin <I> + sin [3 cos $] (26) 
£13 = Xlf sin 8 sin {3 — jlf [cos 8 sin (3 cos $ + cos (3 sin $] 

— Zlf [cos 8 sin/? sin $ — cos/? cos $] (27) 



or, equivalently: 



Xlf = ai cos 8 — 3.2 sin 8 cos (3 + a.3 sin 8 sin (3 (28) 
Ylf = ai sin 8 cos $ + a2 [cos 8 cos f3 cos $ — sin /3 sin <I>] 

—3.3 [cos 8 sin /3 cos <& + cos /3 sin <i>] (29) 
Zlf = ai sin 8 sin <i> + a2 [cos 8 cos /3 sin $ + sin /? cos <&] 

—3.3 [cos 8 sin /3 sin <I> — cos (3 cos $] (30) 



17.2 Orientation of the Incident Beam in the Target Frame 

Under some circumstances, one may wish to specify the target orientation such that xlf (the direction 
of propagation of the radiation) and ylf (usually the first polarization direction) and zlf (= xlf x 
Ylf) refer to certain directions in the TF. Given the definitions of the LF and TF above, this is simply 
an exercise in coordinate transformation. For example, one might wish to have the incident radiation 
propagating along the (1,1,1) direction in the TF (example 14 below). Here we provide some selected 
examples: 

1. XLF = ai, Ylf = £12, zlf = a3 : 8 = 0, $ + /3 = 

2. XLF = ai, YLF = a3, ^lf = -a2 : 8 = 0, $ + /3 = 90° 
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3. 


xlf 


= a2, ylf = 


ai, zlf = 


-as : e = 90°, /3 = 180°, $ = 


4. 


xlf 


= a2, ylf = 


as, ZLF = 


ai : e = 90°, P = 180°, $ = 90° 


5. 


xlf 


= as, YLF = 


ai, ZLF = 


a2 : e = 90°, (3 = -90°, $ = 


6. 


xlf 


= as, YLF = 


a2, ZLF = 


-ai : e = 90°, P = -90°, $ = -90° 


7. 


XLF 


= -ai, YLF 


= a2, ZLF 


= -as : e = 180°, /3 + $ = 180° 


8. 


Xlf 


= -ai, ylf 


= as, ZLF 


= a2 : e = 180°, /3 + $ = 90° 


9. 


Xlf 


= -a2, ylf 


= ai, Zlf 


= as : e = 90°, /3 = 0, $ = 


10. 


Xlf 


= -a2, ylf 


= as, ZLF 


= -ai : e = 90°, /3 = 0, $ = -90° 


11. 


Xlf 


= -as, Ylf 


= ai, ZLF 


= -a2 : e = 90°, P = -90°, $ = 


12. 


Xlf 


= -as, Ylf 


= a2, ZLF 


= ai : e = 90°, P = -90°, $ = 90° 


13. 


Xlf 


= (ai + a2)/V2, YLF = 


= as, ZLF = (ai - a2)/V2 : 9 = 45°, P : 


14. 


Xlf 


= (ai + a2 - 


1- as)/ ^3, YLF = (ai - £12)/ ^2, zlf = (ai + a2 - 



2a3)/V6 : 
e = 54.7356°, P = 135°, $ = 30°. 

17.3 Sampling in 0, $, and P 

The present version, DDSCAT 7.0, chooses the angles /3, 6, and $ to sample the intervals (BETAMI , BETAMX), 
(THETMI , THETMX) , (PHIMIN, PHIMAX), where BETAMI, BETAMX, THETMI, THETMX, PHIMIN, 
PHIMAX are specified in ddscat . par . The prescription for choosing the angles is to: 

• uniformly sample in P; 

• uniformly sample in $; 

• uniformly sample in cos Q. 

This prescription is appropriate for random orientation of the target, within the specified limits of P, $, 
ande. 

Note that when DDSCAT 7.0 chooses angles it handles P and $ differently from Q. The range for P 
is divided into NBETA intervals, and the midpoint of each interval is taken. Thus, if you take BETAMI=0, 
BETAMX=90,NBETA=2 you will get /3 = 22.5° and 67.5°. Similarly, if you take PHIMIN=0, PHIMAX=180, 
NPHI=2 you will get $ = 45° and 135°. 

Sampling in is done quite differently from sampling in P and $. First, as already mentioned above, 
DDSCAT 7.0 samples uniformly in cos Q, not Q. Secondly, the sampling depends on whether NTHETA 
is even or odd. 

• If NTHETA is odd, then the values of 8 selected include the extreme values THETMI and THETMX; 
thus, THETMI=0, THETMX=90, NTHETA=3 will give you 6 = 0, 60°, 90°. 

• If NTHETA is even, then the range of cos 8 will be divided into NTHETA intervals, and the mid- 
point of each interval will be taken; thus, THETMI=0, THETMX=90, NTHETA=2 will give you 
e = 41.41° and 75.52° [cos 8 = 0.25 and 0.75]. 

The reason for this is that if odd NTHETA is specified, then the "integration" over cos 8 is performed 
using Simpson's rule for greater accuracy. If even NTHETA is specified, then the integration over cos 8 
is performed by simply taking the average of the results for the different 8 values. 

If averaging over orientations is desired, it is recommended that the user specify an odd value of 
NTHETA so that Simpson's rule will be employed. 



18 Orientational Averaging 

DDSCAT has been constructed to facilitate the computation of orientational averages. How to go about 
this depends on the distribution of orientations which is applicable. 
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18.1 Randomly-Oriented Targets 

For randomly-oriented targets, we wish to compute the orientational average of a quantity &, $): 

-t pin />! pin 

{Q)^-^J dp J dcosQ J d$0(/3,e,$) . (31) 

To compute such averages, all you need to do is edit the file ddscat .par so that DDSCAT knows 
what ranges of the angles (3, Q, and $ are of interest. For a randomly-oriented target with no symmetry, 
you would need to let f3 run from to 360°, 6 from to 180°, and <!> from to 360°. 

For targets with symmetry, on the other hand, the ranges of (3, 8, and <I> may be reduced. First 
of all, remember that averaging over <i> is relatively "inexpensive", so when in doubt average over to 
360°; most of the computational "cost" is associated with the number of different values of (/3,0) which 
are used. Consider a cube, for example, with axis ai normal to one of the cube faces; for this cube (3 
need run only from to 90°, since the cube has fourfold symmetry for rotations around the axis ai. 
Furthermore, the angle 8 need run only from to 90°, since the orientation (/3,8,$) is indistinguishable 
from (13, 180° - 8, 360° - 

For targets with symmetry, the user is encouraged to test the significance of /3,8,<E> on targets with 
small numbers of dipoles (say, of the order of 100 or so) but having the desired symmetry. 

18.2 Nonrandomly- Oriented Targets 

Some special cases (where the target orientation distribution is uniform for rotations around the x axis = 
direction of propagation of the incident radiation), one may be able to use DDSCAT 7.0 with appropriate 
choices of input parameters. More generally, however, you will need to modify subroutine ORIENT to 
generate a list of NBETA values of /?, NTHETA values of 8, and NPHI values of $, plus two weight- 
ing arrays WGTA ( 1-NTHETA, 1-NPHI ) and WGTB ( 1-NBETA) . Here WGTA gives the weights which 
should be attached to each (8,4>) orientation, and WGTB gives the weight to be attached to each (3 ori- 
entation. Thus each orientation of the target is to be weighted by the factor WGTA x WGTB. For the case 
of random orientations, DDSCAT chooses 8 values which are uniformly spaced in cos 8, and f] and $ 
values which are uniformly spaced, and therefore uses uniform weights 
WGTB=1./NBETA 

When NTHETA is even, DDSCAT sets 

WGTA=1./(NTHETAXNPHI) 

but when NTHETA is odd, DDSCAT uses Simpson's rule when integrating over 8 and 
WGTA= (1/3 or 4/3 or2/3)/(NTHETAxNPHl) 
Note that the program structure of DDSCAT may not be ideally suited for certain highly oriented 
cases. If, for example, the orientation is such that for a given $ value only one 8 value is possible 
(this situation might describe ice needles oriented with the long axis perpendicular to the vertical in the 
Earth's atmosphere, illuminated by the Sun at other than the zenith) then it is foolish to consider all the 
combinations of 8 and <E> which the present version of DDSCAT is set up to do. We hope to improve 
this in a future version of DDSCAT. 

19 Target Generation: Isolated Finite Targets 

DDSCAT contains routines to generate dipole arrays representing finite targets of various geometries, 
including spheres, ellipsoids, rectangular solids, cylinders, hexagonal prisms, tetrahedra, two touching 
ellipsoids, and three touching ellipsoids. The target type is specified by variable CSHAPE on line 9 of 
ddscat .par, up to 12 target shape parameters (SHPARi, SHPAR2, SHPAR3, ...) on line 10. The 
target geometry is most conveniently described in a coordinate system attached to the target which we 
refer to as the "Target Frame" (TF), with orthonormal unit vectors xtf, Ytf, ztf = xtf x yxF- Once 
the target is generated, the orientation of the target in the Lab Frame is accomplished as described in iJlTl 
Every target generation routine will specify 
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• The "occupied" lattice sites; 

• The composition associated with each occupied lattice site; 

• Two "target axes" ai and a2 that are used as references when specifying the target orientation; and 

• The location of the Target Frame origin of coordinates. 
Target geometries currently supported include: 

19.1 FROM_FILE = Target composed of possibly anisotropic material, defined 
by list of dipole locations and "compositions" obtained from a file 

If anisotropic, the "microcrystals" in the target are assumed to be aligned with the principal axes of 
the dielectric tensor parallel to xtf, Ytf, and ztf- This option causes DDSCAT to read the target 
geometry and composition information from a file shape . dat instead of automatically generating 
one of the geometries for which DDSCAT has built-in target generation capability. The shape . dat 
file is read by routine REASHP (file reashp . f 90). The file shape . dat gives the number N of 
dipoles in the target, the components of the "target axes" ai and 3.2 in the Target Frame (TF), the vector 
a;o(l — 3) determining the correspondence between the integers IXYZ and actual coordinates in the TF, 
and specifications for the location and "composition" of each dipole. The user can customize REASHP as 
needed to conform to the manner in which the target description is stored in file shape . dat. However, 
as supplied, REASHP expects the file shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements 

• N = number of dipoles in target 

• o,ix fliy o-iz = x,y,z components (in TF) of ai 

• CL2x a2y o,2z = x,y,z Components (in TF) of a2 

• dx/d dy/d dz/d = 1. 1. 1. = relative spacing of dipoles in xtf, Ytf, ztf directions 

• xqx xoy xq^ = TF coordinates xtf/c^ UTp/d ztf/c? corresponding to lattice site IXYZ= 

• (line containing comments) 

• dummy IXYZ (1, 1) IXYZ (1,2) IXYZ (1,3) IC0MP(1,1) IC0MP(1,2) IC0MP(1,3) 

• dummy IXYZ (2, 1) IXYZ (2, 2) IXYZ (2, 3) IC0MP(2,1) ICOMP(2,2) ICOMP(2,3) 

• dummy IXYZ (3, 1) IXYZ (3, 2) IXYZ (3, 3) IC0MP(3,1) ICOMP(3,2) ICOMP(3,3) 

• ... 

• dummy IXYZ (J, 1) IXYZ (J, 2) IXYZ (J, 3) ICOMP (J, 1) ICOMP(J, 2) IC0MP(J,3) 

• ... 

• dummy IXYZ (N, 1 ) IXYZ (N, 2) IXYZ (N, 3) ICOMP (N,l) ICOMP (N, 2) ICOMP (N, 3) 

where dummy is a string that does not contain spaces, commas, or other separators. 

If the target material at location J is isotropic, ICOMP ( J, 1 ) , ICOMP ( J, 2 ) , and ICOMP ( J, 3 ) have 

the same value. 

19.2 ANIFRMFIL = General anistropic target defined by list of dipole locations, 
"compositions", and material orientations obtained from a file 

This option causes DDSCAT to read the target geometry information from a file shape . dat instead 
of automatically generating one of the geometries listed below. The file shape . dat gives the num- 
ber of dipoles in the target, the components of the "target axes" ai and 3.2 in the Target Frame 
(TF), the vector a::o(l — 3) determining the correspondence between the integers IXYZ and actual co- 
ordinates in the TF, and specifications for the location and "composition" of each dipole. For each 
dipole J, the file shape . dat provides the location IXYZ (J, 1-3 ) , the composition identifier integer 
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ICOMP (J, 1-3 ) specifying the dielectric function corresponding to the three principal axes of the di- 
electric tensor, and angles 9df, 'J'df, and /3df specifying the orientation of the local "Dielectric Frame" 
(DF) relative to the "Target Frame" (TF) (see ij26] l. The DF is the reference frame in which the dielectric 
tensor is diagonalized. The Target Frame is the reference frame in which we specify the dipole locations. 

The shape . dat file is read by routine REASHP (file reashp . f 90). The user can customize 
REASHP as needed to conform to the manner in which the target geometry is stored in file shape . dat. 
However, as supplied, REASHP expects the file shape . dat to have the following structure: 

• one line containing a description; the first 67 characters will be read and printed in various output 
statements. 

• N = number of dipoles in target 

• o-ix o-iy a-iz = x,y,z components (in Target Frame) of ai 

• o.2x o.2y a2z = x,y,z components (in Target Frame) of a2 

• dx/d dy/d dz/d = 1. 1. 1. = relative spacing of dipoles in xtf, Ytf, ztf directions 

• xqx XQy xqz = TF coordinates .ttf / d utf / d ztf / d corresponding to lattice site I XY Z= 

• (line containing comments) 

• dummy IXYZ {1, 1-3) ICOMP (1,1-3) THETADF(l) PHIDF(l) BETADF(l) 

• dummy IXYZ {2, 1-3) ICOMP (2, 1-3) THETADF(2) PHIDF(2) BETADF(2) 

• dummy IXYZ {3, 1-3) ICOMP (3, 1-3) THETADF(3) PHIDF(3) BETADF(3) 

• ... 

• dummy IXYZ {J, 1-3) ICOMP (J, 1-3) THETADF ( J) PHIDF(J) BETADF(J) 

• ... 

• dummy IXYZ {N, 1-3) ICOMP (N, 1-3) THETADF (N) PHIDF (N) BETADF (N) 

Where dummy is a string that does not contain spaces, commas, or other separators. THETADF PHIDF 
BETADF should be given in degrees. 

19.3 ANIELLIPS = Homogeneous, anisotropic ellipsoid. 

SHPARi, SHPAR2, SHPAR3 define the ellipsoidal boundary; 

f Ztf/^V , f yTF/rf V , f ZTF/d\'^ _1 ^^2-) 



V SHPARi / V SHPAR2 / V SHPAR3 / 4 
The TF origin is located at the centroid of the ellipsoid. 

19.4 ANI_ELL_2 = Two touching, homogeneous, anisotropic ellipsoids, with dis- 
tinct compositions 

Geometry as for ELLIPS0_2; SHPARi, SHPAR2, SHPAR3 have same meanings as for ELLIPS0_2. 

Target axes ai = (1, 0, 0)tf and 3.2 = (0, 1, 0)tf- 

Line connecting ellipsoid centers is 1 1 ai = xtf ■ 

TF origin is located between ellipsoids, at point of contact. 

It is assumed that (for both ellipsoids) the dielectric tensor is diagonal in the TF. User must set NC0MP=6 
and provide xx, yy, zz components of dielectric tensor for first ellipsoid, and xx, yy, zz components of 
dielectric tensor for second ellipsoid (ellipsoids are in order of increasing xtf)- 
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19.5 ANI_ELL_3 = Three touching homogeneous, anisotropic ellipsoids with 
same size and orientation but distinct dielectric tensors 

SHPARi, SHPAR2, SHPAR3 have same meanings as for ELLIPS0_3. 

Target axis ai = (1, 0, 0)tf (along line of ellipsoid centers), and 3.2 ~ (0, 1, 0)tf- 

TF origin is located at center of middle ellipsoid. 

It is assumed that dielectric tensors are all diagonal in the TF. User must set NC0MP=9 and provide xx, 
yy, zz elements of dielectric tensor for first ellipsoid, xx, yy, zz elements for second ellipsoid, and xx, 
yy, zz elements for third ellipsoid (ellipsoids are in order of increasing x'tf)- 

19.6 ANIRCTNGL = Homogeneous, anisotropic, rectangular solid 

X, y, z lengths/rf = SHPARi, SHPAR2, SHPAR3. 

Target axes ai = (1, 0, 0)tf and 3.2 = (0, 1, 0)tf in the TF. 

(iTF, 2/TF, ^tf) = (0, 0, 0) at middle of upper target surface, (where "up" = xtf)- (The target surface 

is taken to be d/2 about the upper dipole layer.) 

Dielectric tensor is assumed to be diagonal in the target frame. 

User must set NC0MP=3 and supply names of three files for e as a function of wavelength or energy: first 
for €xx, second for Eyy, and third for e^z, as in the following sample ddscat . par file: 

' =============-^-^^^ Parameter file ' 

PRELIMINARIES ****' 

'NOTORQ'= CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGS2'= CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT'= CMETHD*6 (GPFAFT, FFTMKL) 
'GKDLDR'= CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN'= CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 
' ANIRCTNGL' = CSHAPE*9 shape directive 

Initial Memory Allocation 
10 25 50 = upper bound on target extent 

Target Geometry and Composition 
'ANIRCTNGL' = CSHAPE*9 shape directive 

10 25 50 = shape parameters SHPARI, SHPAR2, SHPAR3 
3 = NCOMP = number of dielectric materials 

' /u/draine/DDA/diel/m2 . 0_0 . 1 ' = name of file containing dielectric function 
' /u/draine/DDA/diel/ml . 5 0_0 .00' 
' /u/draine/DDA/diel/ml . 5 0_0 .00' 
'++++ Error Tolerance ****' 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

' **** Interaction cutoff parameter for PBC calculations ****' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

' **** Angular resolution for calculation of <cos>, etc. ****' 

2.0 = ETASCA (number of angles is proportional to [ (2+x) /ETASCA] "2 ) 

' **** Wavelengths (micron) ***★' 

500. 500. 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
' **** Effective Radii (micron) **** ' 

30.60 30.60 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' **** Define Incident Polarizations 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

' **** Specify which output files to write ****' 

1 = IWRKSC (-0 to suppress, -1 to write ".sea" file for each target orient. 

1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA,THETA) 
' Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 
0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 
'**** Specify first IWAV, TRAD, lORI (normally 0) ****' 
= first IWAV, first IRAD, first TORI (0 to begin fresh) 
' **** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

11 12 21 22 31 41 - indices ij of elements to print 
' **** Specify Scattered Directions ****' 

'LFRAME' = CMDFRM*6 ('LFRAME' or 'TFRAME' for Lab Frame or Target Frame) 

2 = NPLANES = number of scattering planes 
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0. 0. 180. 30 = phi, thetan_min, thetan_max, dtheta (in degrees) for plane A 
90. 0. 180. 30 = phi, ... for plane B 



19.7 CONELLIPS = Two concentric ellipsoids 

SHPARi, SHPAR2, SHPAR3 = lengths/d of the oMfer ellipsoid along the xtf, Ytf, ztf axes; 
SHPAR4, SHPAR5, SHPARg = lengths/dof the inner elUpsoid along the xtf, Ytf, ztf axes. 
Target axes ai = (1, 0, 0)tf, sli = (0, 1, 0)tf- 
TF origin is located at centroids of ellipsoids. 

The "core" within the inner ellipsoid is composed of isotropic material 1 ; the "mantle" between inner 
and outer ellipsoids is composed of isotropic material 2. 

User must set NC0MP=2 and provide dielectric functions for "core" and "mantle" materials. 



19.8 CYLINDERl = Homogene( 

SHPARi = length/d, SHPAR2 = diameter/d, 
SHPAR3 = 1 for cylinder axis ai || xtf: ai 
SHPAR3 = 2 for cylinder axis ai |j Ytf- ai 
SHPAR3 = 3 for cylinder axis ai || ztf^ ai 
TF origin is located at centroid of cylinder. 
User must set NC0MP=1. 



JUS, isotropic finite cylinder 

with 

= (l,0,0)TFanda2 = (0,l,0)TF; 
= (0,l,0)TFand aa = (0,0,1)tf; 
= (0, 0, 1)tf and aa = (1, 0, 0)tf in the TF 



19.9 CYLNDRCAP = Homogeneous, isotropic finite cylinder with hemispheri- 
cal endcaps. 

SHPARi = cylinder length/d (not including end-caps!) and SHPAR2 = cylinder diameter/d, with cylinder 

axis =3.1 = (1,0,0)tf and 3.2 = (0, 1,0)tf- The total length along the target axis (including the 

endcaps) is (SHPARi+SHPAR2)d. 

TF origin is located at centroid of cylinder. 

User must set NC0MP=1. 



19.10 DSKRCTNGL = Disk on top of a homogeneous rectangular slab 

This option causes DDSCAT to create a target consisting of a disk of composition 1 resting on top of a 
rectangular block of composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic, 
ddscat .par should set NCOMP to 2 . 

The cylindrical disk has thickness SHPARi xd in the x-direction, and diameter SHPARaxd. The rect- 
angular block is assumed to have thickness SHPAR3xd in the x-direction, length SHPAR4X(i in the 
y-direction, and length SHPARsxd in the z-direction. The lower surface of the cylindrical disk is in the 
a; = plane. The upper surface of the slab is also in the x — plane. 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the a; = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at {x/d, y/d, z/d) ~ {j^ + 0.5, jy + ^y,jz + A^), where 
jx, jy, and jz are integers. Ay = or 0.5 depending on whether SHPAR4 is even or odd. = or 0.5 
depending on whether SHPAR5 is even or odd. 

Dipoles representing the disk are located at 
x/d = 0.5, 1.5, [int(SHPAR4 + 0.5) - 0.5] 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
Oeff = (3VTuc/47r)^/'^, where Vtuc is the total volume of solid material in one TUC. For this geom- 
etry, the number of dipoles in the target will be approximately N = [SHPARi x SHPAR2 x SHPAR3 + 
(7r/4)((SHPAR4)-^ x SHPAR5)], although the exact number may differ because of the dipoles are required 
to be located on a rectangular lattice. The dipole spacing d in physical units is determined from the 
specified value of and the number N of dipoles in the target: d = (47r/3iV)^/'^acff . This option 
requires 5 shape parameters: 
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The pertinent line in dd scat. par should read 

SHPARl SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARl = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR2 = (disk diameter)/d 

SHPAR3 = (brick thickness in xtf directionVd [material 2] 
SHPAR4 = (brick thickness in yxF directionVd 
SHPAR5 = (brick thickness in ztf direction)/d 



The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters ( S HP ARl + S HP AR4 ) , 
SHPAR2, and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
and SHPAR5. 

The physical size of the TUC is specified by the value of a^s (in physical units, e.g. cm), specified in the 
file ddscat . par. 

The "computational volume" is determined by (SHPARi + SHPAR4)x SHPAR2X SHPAR3. 

The target axes (in the TF) are set to ai = xtf = (I7 0, 0)tf - i-e-. normal to the "slab" - and 
^2 = Ytf = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all other 
targets, set by the usual orientation angles /?, Q, and $ (see |T7] above); for example, = would be 
for radiation incident normal to the slab. 



19.11 DW1996TAR = 13 block target used by Draine & Weingartner (1996). 

Single, isotropic material. Target geometry was used in study by Draine & Weingartner (1996) of radia- 
tive torques on irregular grains, ai and a2 are principal axes with largest and second-largest moments of 
inertia. User must set NC0MP = 1. Target size is controlled by shape parameter SHPAR ( 1 ) = width of 
one block in lattice units. 
TF origin is located at centroid of target. 



19.12 ELLIPSOID = Homogeneous, isotropic ellipsoid. 

"Lengths" SHPARi, SHPAR2, SHPAR3 in the x, y, z directions in the TF: 

/ .TTF \ ^ ^ / j/TF V ^ / ^:tf 
V SHPARl VSHPAR2d/ VSHPAR3(i 

where d is the interdipole spacing. 

The target axes are set to ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf- 
Target Frame origin = centroid of ellipsoid. 
User must set NC0MP=1 on line 9 of ddscat . par. 
A homogeneous, isotropic sphere is obtained by setting SHPARi = SHPAR2 = SHPAR3 = diameter/d. 




19.13 ELLIPSO_2 = Two touching, homogeneous, isotropic ellipsoids, with dis- 
tinct compositions 

SHPARl, SHPAR2, SHPAR3=x-length/(i, y-length/d, z-length/d of one ellipsoid. The two ellipsoids 
have identical shape, size, and orientation, but distinct dielectric functions. The line connecting ellipsoid 
centers is along the XxF-axis. Target axes ai ~ (1,0,0)tf [along line connecting ellipsoids] and 
a2 = (0,1,0)tf. 

Target Frame origin = midpoint between ellipsoids (where ellipsoids touch). 

User must set NC0MP=2 and provide dielectric function file names for both ellipsoids. Ellipsoids are 
in order of increasing xtf: first dielectric function is for ellipsoid with center at negative xtf, second 
dielectric function for elUpsoid with center at positive xtf • 
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19.14 ELLIPSO_3 = Three touching homogeneous, isotropic ellipsoids of equal 
size and orientation, but distinct compositions 

SHPARi, SHPAR2, SHPAR3 have same meaning as forELLIPS0_2. Line connecting ellipsoid centers is 
parallel to xtf axis. Target axis ai = (1, 0, 0)tf (along line of ellipsoid centers), and 3.2 = (0, 1, 0)tf- 
Target Frame origin = centroid of middle ellipsoid. 

User must set NC0MP=3 and provide (isotropic) dielectric functions for first, second, and third ellipsoid. 



19.15 HEX_PRISM = Homogeneous, isotropic hexagonal prism 

SHPARi = (Length of prism)/rf = (distance between hexagonal faces)/(i, 

SHPAR2 = (distance between opposite vertices of one hexagonal face)/(i = 2xhexagon side/d. 

SHPAR3 selects one of 6 orientations of the prism in the Target Frame (TF). 

Target axis ai is along the prism axis (i.e., normal to the hexagonal faces), and target axis 3.2 is normal 
to one of the rectangular faces. There are 6 options for SHPAR3: 

SHPAR3 = 1 for ai II Xtf and 3.2 || yxF ; SHPAR3 = 2 for ai || xtf and a2 || ztf ; 
SHPAR3 = 3 for ai II yxF and a.2 \\ xtf ; SHPAR3 = 4 for ai || yxF and a2 || ztf ; 
SHPAR3 = 5 for ai II Ztf and a.2 \\ xtf ; SHPAR3 = 6 for ai || ztf and a2 || yxF 
TF origin is located at the centroid of the target. 
User must set NCOMP=L 

19.16 LAYRDSLAB = Multilayer rectangular slab 

Multilayer rectangular slab with overall x, y, z lengths Ux — SHPARi x d 

Qy = SHPAR2 X rf, 

= SHPAR3 X d. 

Upper surface is at xtf = 0, lower surface at xtf = —SHPARi x d 

SHPAR4 = fraction which is composition 1 (top layer). 

SHPAR5 = fraction which is composition 2 (layer below top) 

SHPARg = fraction which is composition 3 (layer below comp 2) 

1 - (SHPAR4 + SHPAR5 + SHPARe) = fraction which is composition 4 (bottom layer). 

To create a bilayer slab, just set SHPAR5 = SHPARe = 

To create a trilayer slab, just set SHPARe = 

User must set NCOMP=2,3, or 4 and provide dielectric function files for each of the two layers. Top 
dipole layer is at xtf = —d/2. Origin of TF is at center of top surface. 

19.17 MLTBLOCKS = Homogeneous target constructed from cubic "blocks" 

Number and location of blocks are specified in separate file blocks . par with following structure: 
one line of comments (may be blank) 
PRIN (= or 1 - see below) 
N (= number of blocks) 
B (= width/rf of one block) 

Xtf Utf Ztf (= position of 1st block in units of Be?) 
Xtf Utf ztf (= position of 2nd block in units of Bd) ) 

Xtf Utf ztf (= position of Nth block in units of Be?) 
If PRIN=0, then ai = (1, 0, 0)tf, ^2 = (0, 1, 0)tf- If PRIN=1, then ai and £12 are set to principal axes 
with largest and second largest moments of inertia, assuming target to be of uniform density. User must 

set NC0MP=1. 
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19.18 RCTGLPRSM = Homogeneous, isotropic, rectangular solid 

X, y, z lengths/d = SHPARi, SHPAR2, SHPAR3. 
Target axes ai = (1,0, 0)tf and a2 = (0, 1, 0)tf- 

TF origin at center of upper surface of solid: target extends from xtf /d = — SHPARi to 0, 
j/Tp/d from -0.5 x SHPAR2 to +0.5 x SHPAR2 
ZTp/d from -0.5 x SHPAR3 to +0.5 x SHPAR3 
User must set NC0MP=1. 

19.19 RCTGLBLK3 = Stack of 3 rectangular blocks, with centers on the xtf 
axis. 

Each block consists of a distinct material. There are 9 shape parameters: 

SHPARi = (upper solid thickness in xtf direction)/^ [material 1] 

SHPAR2 = (upper solid width in yxp direction)/^ 

SHPAR3 = (upper solid width in ztf direction)/^ 

SHPAR4 = (middle solid thickness in xtf direction)/d [material 2] 

SHPAR5 = (middle solid width in ytf direction)/d 

SHPARg = (middle solid width in ztf direction)/(i 

SHPAR7 = (lower solid thickness in xtf direction)/^ [material 3] 

SHPARg = (lower solid width in yxF direction)/^ 

SHPARg = (lower solid width in ztf direction)/^ 

TF origin is at center of top surface of material 1 . 

19.20 SPHERES_N = Multisphere target = union of N spheres of single isotropic 
material 

Spheres may overlap if desired. The relative locations and sizes of these spheres are defined in an 
external file, whose name (enclosed in single quotes) is passed through ddscat .par. The length of 
the file name should not exceed 80 characters. The pertinent line in ddscat . par should read 
SHPARi SHPAR2 '_^/ename ' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR2= to have oi = (1, 0, 0)tf, 02 = (0, 1, 0)tf- 

SHPAR2= 1 to use principal axes of moment of inertia tensor for ai (largest /) and 02 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPAR ( 1 ) , which is the extent of the multisphere target in the x-direction, in units of the lattice spacing 
d. The file filename ' should have the following structure: 

N (= number of spheres) 

one line of comments (may be blank) 

Xi yi z\ fli (arb. units) 

X2 j/2 22 a2 (arb. units) 

a^iv Vn zn Oat (arb. units) 
where Xj, yj, zj are the coordinates (in the TF) of the center, and aj is the radius of sphere j. 
Note that xj, yj, Zj, aj (j = 1, ■■■,N) establish only the shape of the A^— sphere target. For instance, 
a target consisting of two touching spheres with the line between centers parallel to the x axis could 
equally well be described by lines 3 and 4 being 

00 0.5 

1 00 0.5 

or 

000 1 
200 1 

The actual size (in physical units) is set by the value of Ucs specified in dds cat . par, where, as always. 



TARGET GENERATION: ISOLATED FINITE TARGETS 



Ooff = {SV/AttY^^, where V is the total volume of material in the target. 

Target axes ai and a.2 are set to be principal axes of moment of inertia tensor (for uniform density), 
where ai corresponds to the largest eigenvalue, and 3.2 to the intermediate eigenvalue. 
The TF origin is taken to be located at the volume-weighted centroid. 
User must set NC0MP=1. 

19.21 SPHROID_2 = Two touching homogeneous, isotropic spheroids, with dis- 
tinct compositions 

First spheroid has length SHPARi along symmetry axis, diameter SHPAR2 perpendicular to symmetry 
axis. Second spheroid has length SHPAR3 along symmetry axis, diameter SHPAR4 perpendicular to 
symmetry axis. Contact point is on line connecting centroids. Line connecting centroids is in xtf 
direction. Symmetry axis of first spheroid is in yxF direction. Symmetry axis of second spheroid is in 
direction yxF cos(SHPAR5) + ztf sin(SHPAR5), and SHPAR5 is in degrees. If SHPARg = 0., then target 
axes ai = (1, 0, 0)tf, 3-2 = (0, 1, 0)tf- If SHPARg = 1., then axes ai and 3.2 are set to principal axes 
with largest and 2nd largest moments of inertia assuming spheroids to be of uniform density. 
Origin of TF is located between spheroids, at point of contact. 
User must set NC0MP=2 and provide dielectric function files for each spheroid. 

19.22 SPH_ANI_N = Multisphere target consisting of the union of spheres of 
various materials, possibly anisotropic 

Spheres may NOT overlap. The relative locations and sizes of these spheres are defined in an external 
file, whose name (enclosed in single quotes) is passed through ddscat .par. The length of the file 
name should not exceed 80 characters. Target axes ai and 3.2 are set to be principal axes of moment 
of inertia tensor (for uniform density), where ai corresponds to the largest eigenvalue, and a.2 to the 
intermediate eigenvalue. 

The TF origin is taken to be located at the volume-weighted centroid. 
The pertinent line in ddscat . par should read 
SHPARi SHPAR2 '^/ename ' (quotes must be used) 

where SHPARi = target diameter in x direction (in Target Frame) in units of d 
SHPAR2= to have ai = (1, 0, 0)tf, 02 = (0, 1, 0)tf in Target Frame. 

SHPAR2= 1 to use principal axes of moment of inertia tensor for ai (largest /) and a2 (intermediate /). 
filename is the name of the file specifying the locations and relative sizes of the spheres. 
The overall size of the multisphere target (in terms of numbers of dipoles) is determined by parameter 
SHPARi, which is the extent of the multisphere target in the x-direction, in units of the lattice spacing d. 
The file filename ' should have the following structure: 

N (= number of spheres) 

one line of comments (may be blank) 

xi yi zi ri Cxi Cyi Czi Qdf,! *DFa /3df,i 

X2 2/2 Z2 r2 Cx2 Cy2 Cz2 6df,2 <&df,2 /3df,2 

XN yN ZN TN CXN CyN Czn QdF.N <&DF,Af PdF,N 

where xj, yj, Zj are the coordinates of the center, and rj is the radius of sphere j (arbitrary units), 
Cxj, Cyj, Czj are integers specifying the "composition" of sphere j in the x,y,z directions in the 
"Dielectric Frame" (see S|26]l of sphere j, and Qbfj 'i'DF.j /3dfj are angles (in radians) specifying 
orientation of the dielectric frame (DF) of sphere j relative to the Target Frame. Note that xj, yj, Zj, rj 
(j = 1, A^) establish only the shape of the A^— sphere target, just as for target option NSPHER. The 
actual size (in physical units) is set by the value of Ocff specified in ddscat . par, where, as always, 
Ooff = (SV/AttY^^, where V is the volume of material in the target. 

User must set NCOMP to the number of different dielectric functions being invoked (i.e., the range of 

{Cxj,Cyj,Czj}. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file filename ' would be 
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simply 

TV (= number of spheres) 

one line of comments (may be blank) 

a;i yi zi ri 1 1 1 

a;2 ^2 2:2 r2 111 

xn Vn zn 1 1 1000 



19.23 TETRAHDRN = Homogeneous, isotropic tetrahedron 

SHPARi=length/(i of one edge. Orientation: one face parallel to yxF, ztf plane, opposite "vertex" is 
in +xtf direction, and one edge is parallel to zxp. Target axes ai = (1, 0, 0)tf [emerging from one 
vertex] and a2 = (0, 1, 0)tf [emerging from an edge] in the TF. User must set NC0MP=1. 

19.24 TRNGLPRSM = Triangular prism of homogeneous, isotropic material 

SHPARi, SHPAR2, SHPAR3, SHPAR4 = a/d, 6/a, c/a, L/a 

The triangular cross section has sides of width a, b, c. L is the length of the prism, d is the lattice 
spacing. The triangular cross-section has interior angles a, (3, 7 (opposite sides a, h, c) given by cos a = 
(6^ + - a^)/2hc, cos/3 = (a^ + - h'^)/2ac, cosj = (a^ + 6^ - c^)/2ah. In the Target Frame, 
the prism axis is in the x direction, the normal to the rectangular face of width a is (0,1,0), the normal to 
the rectangular face of width h is (0, — cos 7, sin 7), and the normal to the rectangular face of width c is 
(0, — C0S7, — sin7). 

19.25 UNIAXICYL = Homogeneous finite cylinder with uniaxial anisotropic di- 
electric tensor 

SHPARi, SHPAR2 have same meaning as for CYLINDERl. Cylinder axis = ai = (1,0, 0)tf, a2 = 
(0, 1, 0)tf- It is assumed that the dielectric tensor e is diagonal in the TF, with €yy = e^z- User must set 
NC0MP=2. Dielectric function 1 is for E || ai (cylinder axis), dielectric function 2 is for E _L ai. 

19.26 IModifying Existing Routines or Writing New Ones 

The user should be able to easily modify these routines, or write new routines, to generate targets with 
other geometries. The user should first examine the routine target . f 90 and modify it to call any 
new target generation routines desired. Alternatively, targets may be generated separately, and the target 
description (locations of dipoles and "composition" corresponding to x,y,z dielectric properties at each 
dipole site) read in from a file by invoking the option FROM_FILE in ddscat . f 90. 

Note that it will also be necessary to modify the routine reapar . f 9 so that it will accept whatever 
new target option is added to the target generation code . 

19.27 Testing Target Generation using CALLTARGET 

It is often desirable to be able to run the target generation routines without running the entire DDSCAT 
code. We have therefore provided a program CALLTARGET which allows the user to generate targets 
interactively; to create this executable just typ^H 

Non-Linux sites: The source code for CALLTARGET is in the file CALLTARGET . f 90. You must compile and link 

CALLTARGET. f 90, ddcommon . f 90, dsyevj3.f90, errmsg.f90, gasdev.f90, p_lm.f90, prinaxis . f 90, 
ran3.f90, reashp.f90, sizer.f90, tar2el.f90, tar2sp.f90, tar3el.f90, taranirec . f 90, 
tarblocks . f 90, tarcel.f90, tarcyl.f90, tarcylcap . f 90, tarell.f90, target. f90, targspher . f 90, 
tarhex.f90, tarnas.f90, tarnsp.f90, tarpbxn.f90, tarprsm. f 90, tarrctblkS . f 90, tarrecrec . f 90, 
tarslblin . f 90, tartet . f 90, and wrimsg. f 90. 
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make calltarget . 

The program calltarget is to be run interactively; the prompts are self-explanatory. You may need 
to edit the code to change the device number IDVOUT as for DDSCAT (see i]6.7l above). 

After running, calltarget will leave behind an ASCII file target . out which is a list of the 
occupied lattice sites in the last target generated. The format of target . out is the same as the format 
of the shape . dat files read if option FROM_FILE is used (see above). Therefore you can simply 

mv target. out shape.dat 
and then use DDSCAT with the option FROM.FILE (or option ANIFRMFIL in the case of anisotropic 
target materials with arbitrary orientation relative to the Target Frame) in order to input a target shape 
generated by CALLTARGET. 

Note that CALLTARGET - designed to generate finite targets - can be used with some of the "PBC" 
target options (see i|20l below) to generate a list of dipoles in the TUC. At the moment, CALLTARGET 
has support for tai-get options BISLINPBC, DSKBLYPBC, and DSKRCTPBC. 



20 Target Generation: Periodic Targets 

A periodic target consists of a "Target Unit Cell" (TUC) which is then repeated in either the yxF direc- 
tion, the ztf direction, or both. Please see Draine & Flatau (2008) for illustration of how periodic targets 
are assembled out of TUCs, and how the scattering from these targets is constrained by the periodicity. 
The following options are included in DDSCAT: 



20.1 BISLINPBC = Bi-Layer Slab with Parallel Lines 

The target consists of a bi-layer slab, on top of which there is a "line" with rectangular cross-section. 
The "line" on top is composed of material 1, has height Xi (in the xtf direction), width Yi (in the ytf 
direction), and is infinite in extent in the ztf direction. 

The bilayer slab has width Y2 (in the yxF direction). It is consists of a layer of thickness X2 of material 

2, on top of a layer of material 3 with thickness X3. 

SHPARi = Xi /d {Xi = thickness of fine) 

SHPAR2 = Yi /d {Yi = width of line) 

SHPAR3 = X2 1 d {X2 = thickness of upper layer of slab) 

SHPAR4 = X3 /d (X3 = thickness of lower layer of slab) 

SHPAR5 = y2/rf (12 = width of slab) 

SHPARg = Py/d [Py = periodicity in yxF direction). 

If SHPARg = 0, the target is NOT periodic in the yxF direction, consisting of a single column, infinite 
in the ztf direction. 



20.2 CYLNDRPBC = Target consisting of homogeneous cylinder repeated in tar- 
get y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an infinite array of cylinders. The individual 
cylinders are assumed to be homogeneous and isotropic, just as for option RCTNGL (see i^l9.18l l. 
Let us refer to a single cylinder as the Target Unit Cell (TUC). The TUC is then repeated in the target y- 
and/or z-directions, with periodicities PYDxrf and PZDxd, where d is the lattice spacing. To repeat in 
only one direction, set either P YD or P ZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 



SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where SHPARi, SHPAR2, SHPAR3, SHPAR4, SHPAR5 are numbers: 
SHPARi = cylinder length along axis (in units of d) in units of d 
SHPAR2 = cylinder diameter/d 
SHPAR3 = 1 for cylinder axis || xtf 
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= 2 for cylinder axis || yxF 

= 3 for cylinder axis || ztf (see below) 
SHPAR4 = PYD = periodicity/d in yxF direction ( = to suppress repetition) 
SHPAR5 = PZD = periodicity/d in ztf direction ( = to suppress repetition) 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARl 
and SHPAR2. The orientation of a single cylinder is determined by SHPAR3. The periodicity in the TF 
y and z directions is determined by parameters SHPAR4 and SHPAR5. 

The physical size of the TUC is specified by the value of a^tf (in physical units, e.g. cm), specified in 
the file ddscat . par, with the usual correspondence d = (47r/3iV)^/'^acff, where N is the number of 
dipoles in the TUC. 

With target option CYLNDRPBC, the target becomes a periodic structure, of infinite extent. 

• If NPY > and NPZ > 0, then the target cylindrical TUC repeats in the yxF direction, with 
periodicity NPY x d. 

• If NPY = and NPZ > then the target cylindrical TUC repeats in the ztf direction, with 
periodicity NPZ x d. 

• If NPY > and NPZ > then the target cyHndrical TUC repeats in the yxF direction, with 
periodicity NPY x d, and in the ztf direction, with periodicity NPZ x d. 

Target Orientation: The target axes (in the TF) are set to ai = (1, 0, 0)tf and a2 = (0, 1, 0)tf- Note 
that ai does not necessarily coincide with the cylinder axis: individual cylinders may have any of 3 
different orientations in the TF. 



Example 1: One could construct a single infinite cylinder with the following two lines in ddscat .par; 
100 1 100 

1.0 100.49 2 1.0 0. 

The first line ensures that there will be enough memory allocated to generate the target. The TUC 
would be a thin circular "slice" containing just one layer of dipoles. The diameter of the circular slice 
would be about 100.49rf in extent, so the TUC would have approximately {tt/A) x (100.49)^ = 7931 
dipoles (7932 in the actual realization) within a 100 x 1 x 100 "extended target volume". The TUC would 
be oriented with the cylinder axis in the yxF direction (SHPAR3=2) and the structure would repeat in 
the yxF direction with a period of 1.0 x d. SHPAR5=0 means that there will be no repetition in the z 
direction. As noted above, the "target axis" vector ai = xtf- 

Note that SHPARi, SHPAR2, SHPAR4, and SHPAR5 need not be integers. However, SHPAR3, deter- 
mining the orientation of the cylinders in the TF, can only take on the values 1,2,3. 

The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles (3, 8, and $ (see i |T7] above): for example, 6 = would be for radiation incident 
normal to the periodic structure. 

If either PYD = or PZD ~ (but not both), then the scattering directions are specified as discussed 
in Sj2L2l by specifying a diffraction order AI and the azimuthal angles ^ for which scattering is to be 
calculated for the selected value of M. In many cases M = is the only allowed diffraction order 

If PYD > and PZD > (2-d array of cylinders) then the scattering directions are specified as 
discussed in i l21.3l by specifying the diffraction order {M, N). In many cases, {AI, N) = (0, 0) is the 
only allowed diffraction order For each selected diffraction order one transmitted wave direction and 
one reflected wave direction will be calculated, with the dimensionless 4x4 scattering matrix 5*^2^) 
calculated for each scattering direction. At large distances from the infinite slab, the scattered Stokes 
vector in the (M,N) diffraction order is 
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where /;„ j is the incident Stokes vector. 

20.3 DSKBLYPBC = Target consisting of a periodic array of disks on top of a 
two-layer rectangular slabs. 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of Target Unit 
Cells (TUCs), each TUC consisting of a disk of composition 1 resting on top of a rectangular block con- 
sisting of two layers: composition 2 on top and composition 3 below. Materials 1, 2, and 3 are assumed 
to be homogeneous and isotropic. 
This option requires 8 shape parameters: 
The pertinent line in dd scat. par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARe SHPAR7 SHPARs 
where 

SHPARi = disk thickness in x direction (in Target Frame) in units of d 

SHPAR2 = (disk diameter)/rf 

SHPAR3 = (upper slab thickness)/(i 

SHPAR4 = (lower slab thickness)/^ 

SHPAR5 = (slab extent in yxF direction)/^ 

SHPARg = (slab extent in z^f direction)/^ 

SHPAR7 = period in yxF direction/d 

SHPARg = period in Ztf direction/d 

The physical size of the TUC is specified by the value of a^s (in physical units, e.g. cm), specified 
in the file ddscat . par. 

The "computational volume" is determined by 

(SHPARi + SHPAR3 + SHPAR4)xSHPAR5XSHPAR6. 

The lower surface of the cylindrical disk is in the a; = plane. The upper surface of the slab is also 
in the a; = plane. It is required that SHPAR2 < min(SHPAR4, SHPAR5). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the x = 
plane (the upper surface of the slab, and the lower surface of the disk). 
In the Target Frame, dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPARi + 0.5) - 0.5] 
and at {y, z) values 

2//d = ±0.5,±1.5,...and 

z/d^ ±0.5, ±1.5, ... satisfying 

(y2 + z'^) < (SHPAR2/2)2d2. 

Dipoles representing the rectangular slab are located at {x/d, y/d, z/d) = {jx + 0.5, jy + Ay, jz + 
Az), where j^, jy, and jz are integers. Ay ~ or 0.5 depending on whether SHPAR5 is even or odd. 
Az = or 0.5 depending on whether SHPARg is even or odd. 

The TUC is repeated in the target y- and z-directions, with periodicities SHPARyxd and SHPARg xd. 

As always, the physical size of the target is fixed by specifying the value of the effective radius 
Oefi = (3VTUc/47r)^/'^, where Vtuc is the total volume of sohd material in one TUC. For this geometry, 
the number of dipoles in the target will be approximately 

N = (7r/4) X SHPARi x (SHPAR2)^ + [SHPAR3 + SHPAR4] x SHPAR5 x SHPARg 

although the exact number may differ because of the dipoles are required to be located on a rectangular 
lattice. The dipole spacing d in physical units is determined from the specified value of Oeff and the 
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number N of dipoles in the target: d ~ (47r/3A^)^/'^aofr. 

The target axes (in the TF) are set to ai = xtf = (1, 0, 0)tf - i-e., normal to the "slab" - and 3.2 = 
Ytf = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all other 
targets, set by the usual orientation angles /3, O, and $ (see 23 above); for example, = would be 
for radiation incident normal to the slab. 

The scattering directions are specified by specifying the diffraction order (M, iV); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix S'^^''^ calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is (Draine & Flatau 2008) 



where j is the incident Stokes vector See Draine (& Flatau (2008) for discussion of reflection and 
transmission coefficients. 

20.4 DSKRCTPBC = Target consisting of homogeneous rectangular brick plus 
a disk, extended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of Target Unit Cells. 
Each Target Unit Cell (TUC) consists of a disk of composition 1 resting on top of a rectangular block of 
composition 2. Materials 1 and 2 are assumed to be homogeneous and isotropic. 

The cylindrical disk has thickness SHPARiXc? in the x-direction, and diameter SHPAR2XC?. The rect- 
angular block is assumed to have thickness SHPARaxd in the x-direction, extent SHPAR4X(i in the 
y-direction, and extent SHPARsxd in the z-direction. The lower surface of the cylindrical disk is 
in the x = plane. The upper surface of the slab is also in the a: = plane. It is required that 
SHPAR2<min(SHPAR4, SHPAR5). 

The Target Frame origin (0,0,0) is located where the symmetry axis of the disk intersects the a; = 
plane (the upper surface of the slab, and the lower surface of the disk). In the Target Frame, dipoles 
representing the rectangular block are located at {xjd, y/d, z/d) = [j^ + 0.5, + Ay,jz + Az), where 
jx, jy, and jz are integers. Ay = or 0.5 depending on whether SHPAR4 is even or odd. Az = or 0.5 
depending on whether SHPAR5 is even or odd. 
= -[int(SHPAR3 + 0.5)], -1. 

jy = -[int(0.5 X SHPAR4 - 0.5) + 1] , ... , = -[int(SHPAR4 + 0.5) - 0.5] 

y/d = -[int(0.5 x SHPAR4 + 0.5) - 0.5], [int(0.5 x SHPAR4 + 0.5) - 0.5] 

z/d ^ -[int(0.5 x SHPAR5 + 0.5) - 0.5], [int(0.5 x SHPAR5 + 0.5) - 0.5] 
where int(x) is the greatest integer less than or equal to x. Dipoles representing the disk are located at 

x/d = 0.5, 1.5, [int(SHPAR4 + 0.5) - 0.5] 
and at {y, z) values 

y/d = ±0.5, ±1.5, ... and 

z/d = ±0.5, ±1.5, ... satisfying 

(y2 + z'^) < (SHPAR5/2)2d2. 
The TUC is repeated in the target y- and z-directions, with periodicities SHPARe x d and SHPAR7 x d. 
As always, the physical size of the target is fixed by specifying the value of the effective radius Gch = 
(3VTUc/47r)^/'^, where Vtuc is the total volume of solid material in one TUC. For this geometry, 
the number of dipoles in the target will be approximately N = [SHPARi x SHPAR2 x SHPAR3 + 
(7r/4)((SHPAR4)-^ x SHPAR5)], although the exact number may differ because of the dipoles are re- 
quired to be located on a rectangular lattice. The dipole spacing d in physical units is determined from 
the specified value of Ucs and the number N of dipoles in the target: d ~ (47r/3Af)^/'^aoff. This option 
requires 7 shape parameters: 
The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARg SHPAR7 
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where 

SHPARi = [disk thickness (in xtf direction)]/^ [material 1] 
SHPAR2 = (disk diameter)/rf 

SHPAR3 = (brick thickness in xtf direction)/^ [material 2] 
SHPAR4 = (brick length in yxF direction)/d 
SHPAR5 = (brick length in ztf direction)/(i 
SHPARg = periodicity in yxF direction/d 
SHPAR7 = periodicity in Ztf direction/d 

The overall extent of the TUC (the "computational volume") is determined by parameters ( SHPARI + SHPAR4 ) , 
max ( SHPAR2 , SHPAR4, and max ( SHPAR3 , SHPAR5 ) . The periodicity in the TF y and z directions is 
determined by parameters SHPARg and SHPAR7. 

The physical size of the TUC is specified by the value of Qcs (in physical units, e.g. cm), specified in the 
file ddscat . par. 

The target is a periodic structure, of infinite extent in the target y- and z- directions. The target axes 
(in the TF) are set to ai = xtf — (1,0, 0)tf - i-S-, normal to the "slab" - and 3.2 = ytf = (0, 1, 0)tf- 
The orientation of the incident radiation relative to the target is, as for all other targets, set by the usual 
orientation angles /?, 0, and <!> (see SjTT] above); for example, 8 = would be for radiation incident 
normal to the slab. 

The scattering directions are specified by specifying the diffraction order {M, N); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix 5'(^'') calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is (Draine & Flatau 2008) 



where lin.j is the incident Stokes vector See Draine & Flatau (2008) for discussion of reflection and 
transmission coefficients. 

20.5 HEXGONPBC = Target consisting of homogeneous hexagonal prism re- 
peated in target y and/or z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic or biperiodic array of hexagonal 
prisms. The individual prisms are assumed to be homogeneous and isotropic, just as for option RCTNGL 
(see 

Let us refer to a single hexagonal prism as the Target Unit Cell (TUC). The TUC is then repeated in 
the target y- and z-directions, with periodicities PYDxd and PZDxd, where d is the lattice spacing. To 
repeat in only one direction, set either P YD or P ZD to zero. 

This option requires 5 shape parameters: The pertinent line in ddscat . par should read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 

whereSHPARi, SHPAR2, SHPAR3, SHPAR4, SHPAR5 are numbers: 
SHPARi = prism length along prism axis (in units of d) in units of d 
SHPAR2 = 2xlength of one hexagonal side/d 

SHPAR3 = 1,2,3,4,5 or 6 to specify prism orientation in the TF (see below) 
SHPAR4 = PYD = periodicity in TF y direction/d 
SHPAR5 = PZD = periodicity in TF z direction/d 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the TF y and z directions is determined by parameters SHPAR4 
and SHPAR5. 

The physical size of the TUC is specified by the value of Ooff (in physical units, e.g. cm), specified in 
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the file ddscat . par, with the usual correspondence d = (47r/3A^)^/'^acff, where N is the number of 
dipoles in the TUC. 

With target option HEXGONPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions (assuming both NP Y and NP Z are nonzero). 

The target axes (in the TF) are set to ai = xtf = (1, 0, 0)tf - i-e., normal to the "slab" - and 
a.2 = yxF = (0, 1,0)tf- 

The individual hexagons may have any of 6 different orientations relative to the slab: Let unit vectors 
h be II to the axis of the hexagonal prism, and let unit vector f be normal to one of the rectangular faces 
of the hexagonal prism. Then 
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= 1 forh 


XTF, f 


yxF 


SHPAR3 
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XTF, f 


Ztf 
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= 6 for h 


XTF, f 
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For example, one could construct a single infinite hexagonal column with the following line in 

ddscat . par: 

2.0 100.0 3 2.0 0. 

The TUC would be a thin hexagonal "slice" containing two layers of dipoles. The edges of the hexagon 
would be about 50c? in extent, so the TUC would have approximately (3v^/2) x 50^ x 2 = 12990 
dipoles (13024 in the actual realization) within a 90 x 2 x 100 "extended target volume". The TUC 
would be oriented with the hexagonal axis in the yxF direction, with zxf normal to a rectangular faces 
of the prism (SHPAR3=3), and the structure would repeat in the yxF direction with a period of 2 x d 
(SHPAR4=2 . 0). SHPAR5=0 means that there will be no repetition in the zxf direction. 

Note that SHPARi, SHPAR2, SHPAR4, and SHPAR5 need not be integers. However, SHPAR5, deter- 
mining the orientation of the prisms in the TF, can only take on the values 1,2,3,4,5,6. 

Important Note: For technical reasons, PYD and PZD must not be smaller than the "extended" 
target extent in the yTF and ztf directions. When the GPFAFT option is used for the 3 -dimensional 
FFT calculations, the extended target volume always has dimensions/d = 2'^3^5'^, where a, b, and c are 
nonnegative integers, with (dimension/o?) > 1). 

The orientation of the incident radiation relative to the target is, as for all other targets, set by the 
usual orientation angles /?, 9, and $ (see i |T7] above): for example, 6 = would be for radiation incident 
normal to the periodic structure. The scattering directions are specified by specifying the diffraction 
order (A/, N); for each diffraction order one transmitted wave direction and one reflected wave direction 
will be calculated, with the dimensionless 4x4 scattering matrix 5'^^'*' calculated for each scattering 
direction. At large distances from the infinite slab, the scattered Stokes vector in the (M,N) diffraction 
order is 

4 

IscaAM, N)^Y1 S■]''^I^n.J (37) 
1=1 

where /,„ j is the incident Stokes vector 

20.6 LYRSLBPBC = Target consisting of layered slab, extended in target y and 
z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of an array of multilayer bricks, layered in 
the Xtf direction. The size of each brick in the yTF and ztf direction is specified. Up to 4 layers are 
allowed. 

The bricks are repeated in the yTF and Ztf direction with a specified periodicity. If Ly = Py and 
Lz = Pz, then the target consists of a continuous multilayer slab. For this case, it is most economical to 

set Ly/d = L:,/d = Py/d = Fjd = 1. 
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The upper surface of the slab is asssume to be located at .ttf = 0. The lower surface of the slab is at 

Xtf = -Lx = - SHPARlxd. 

The multilayer slab geometry is specified with 9 parameters. The pertinent line in ddscat . par should 
read 

SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 SHPARg SHPAR7 SHPARs SHPARg 
where 

SHPARi = L^/d = (brick thickness in Xtf direction)/^ 
SHPAR2 = Ly/d = (brick extent in yxF direction)/d 
SHPAR3 = L./d = (brick extent in ztf direction)/d 
SHPAR4 = fraction of the slab with composition 1 
SHPAR5 = fraction of the slab with composition 2 
SHPARg = fraction of the slab with composition 3 
SHPAR7 = fraction of the slab with composition 4 
SHPARg = Pij/d = (periodicity in yxF direction)/^ 
SHPARg = Pz/d = (periodicity in ztf direction)/d 

For a slab with only one layer, set SHPAR5 = 0, SHPARg = 0, SHPAR7 = 0. 
For a slab with only two layers, set SHPARg = and SHPAR7 = 0. 
For a slab with only three layers, set SHPAR7 = 0. 

The user must set NCOMP equal to the number of nonzero thickness layers. 

The number of dipoles in one TUC is = nint(SHPARi) x nint(SHPAR2) x nint(SHPAR3). 

The fractions SHPAR2, SHPAR3, SHPAR4, and SHPAR5 must sum to 1. The number of dipoles in 
each of the layers will be integers that are close to nint(SHPARi * SHPAR4), nint(SHPARi * SHPAR5), 
nint(SHPARi * SHPARg), nint(SHPARi * SHPAR7), 

The physical size of the TUC is specified by the value of floff (in physical units, e.g. cm), specified 
in the file ddscat . par. Because of the way Ooff is defined {Ana^^/'i = Nd^), it should be set to 



fleff 



With target option LYRSLBPBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i-c, normal to the 
"slab" - and a2 = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles P, O, and $ (see i?T7l above): for example, = would 
be for radiation incident normal to the slab. 

For this option, there are only two allowed scattering directions, corresponding to transmission and 
specular reflection. DDSCAT will calculate both the transmission and reflection coefficients. 
The last two lines in ddscat . par should appear as in the following example ddscat . par file. This 
example is for a slab with two layers: the slab is 26 dipole layers thick: the first layer comprises 76.92% 
of the thickness, the second layer 23.08% of the thickness. The wavelength is 0.532^m, the thickness is 

= (47r/3)i/37v5''^aefi- (47r/3)i/3(26)2/30.009189 = 0.1300/im. 
The upper layer thickness is 0.2308ia; = 0.0300yum 

ddscat . par is set up to calculate a single orientation: in the Lab Frame, the target is rotated through 
an angle Q = 120°, with $ = 0. In this orientation, the incident radiation is propagating in the 
(—0.5, 0.866, 0) direction in the Target Frame, so that is impinging on target layer 2 (Au). 



' =================== Parameter file ===================' 

'**** PRELIMINARIES ****' 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 

'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 

'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 

'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 

'NOTBIN' = CBINFLAG (ALLBIN, ORIBIN, NOTBIN) 
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'**** Initial Memory Allocation ****' 
26 1 1 = upper bounds on size of TUC 
'**** Target Geometry and Composition ****' 
'LYRSLBPBC = CSHAPE*9 shape directive 

26 0.7692 0.2308 0= shape parameters SHPARl, SHPAR2, SHPAR3, SHPAR4, SHPAR5 

2 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/Eagle_2 00 ' = refractive index 1 

' /u/draine/work/DDA/diel/Au_evap' = refractive index 2 

'**** Error Tolerance ****' 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC| E>-ACA | X> ) / (NORM OF AC|E>) 

Interaction cutoff parameter for PBC calculations ****' 
5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 
'**** Angular resolution for calculation of <cos>, etc. ****' 
0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] "2 ) 
'**** Wavelengths (micron) ****' 

0.5320 0.5320 1 ' INV = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

0.009189 0.009189 1 'LIN' = eff. radii ( fir st , last , how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

Specify which output files to write ****' 
1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

120. 120. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, TORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_i j to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 

Specify Scattered Directions ****' 
'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
0. 0. = (M,N) for scattering 

20.7 RCTGL_PBC = Target consisting of homogeneous rectangular brick, ex- 
tended in target y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a biperiodic array of rectangular bricks. The 
bricks are assumed to be homogeneous and isotropic, just as for option RCTNGL (see ijl9.18l l. 
Let us refer to a single rectangular brick as the Target Unit Cell (TUC). The TUC is then repeated 
in the ?/tf- and zxp-directions, with periodicities PYAEFFxaofr and PZAEFFxooff, where Ucs = 
(3VTUc/47r)^/'^, where Vtuc is the total volume of solid material in one TUC. This option requires 
5 shape parameters: 

The pertinent line in ddscat . par should read 
SHPARl SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where 



SHPARl 
SHPAR2 
SHPAR3 
SHPAR4 



(brick thickness)/^ in the a;TF direction 
(brick thickness)/d in the i/tf direction 
(brick thickness)/(i in the ztf direction 
periodicity/d in the i/tf direction 
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SHPAR5 = periodicity/d in the z^p direction 



The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi, 
SHPAR2, and SHPAR3. The periodicity in the j/tf and ztf directions is determined by parameters 

SHPAR4 and SHPAR5. 

The physical size of the TUC is specified by the value of acs (in physical units, e.g. cm), specified in the 
file ddscat . par. 

With target option RCTGL_PBC, the target becomes a periodic structure, of infinite extent in the 
target y- and z- directions. The target axes (in the TF) are set to ai = (1, 0, 0)tf - i-e., normal to the 
"slab" - and 3.2 = (0, 1, 0)tf- The orientation of the incident radiation relative to the target is, as for all 
other targets, set by the usual orientation angles Q, and $ (see i |T7l above): for example, = would 
be for radiation incident normal to the slab. 

The scattering directions are specified by specifying the diffraction order (M, N); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix 5'^''-' calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is 

4 

hcaAM, N)^Y1 4"''-^"-J- (38) 

where Imj is the incident Stokes vector. See Draine & Flatau (2008) for interpretation of the Sij as 
transmission and reflection efficiencies. 



20.8 RECRECPBC = Rectangular solid resting on top of another rectangular 
solid, repeated periodically in target y and z directions using periodic boundary 
conditions 

The TUC consists of a single rectangular "brick", of material 1, resting on top of a second rectangular 
brick, of material 2. The centroids of the two bricks along a line in the xtf direction. The bricks are as- 
sumed to be homogeneous, and materials 1 and 2 are assumed to be isotropic. The TUC is then repeated 
in the ^tf" and ztf -directions, with periodicities SHPAR4 x d and SHPAR5 x d. is the total volume of 
solid material in one TUC. This option requires 5 shape parameters: 
The pertinent line in ddscat . par should read 



SHPARi SHPAR2 SHPAR3 SHPAR4 SHPAR5 



where 

SHPARi = (upper brick thickness)/^ in the xtf direction 
SHPAR2 = (upper brick thickness)/(i in the yxF direction 
SHPAR3 = (upper brick thickness)/(i in the ztf direction 
SHPAR4 = (lower brick thickness)/(i in the xtf direction 
SHPAR5 = (lower brick thickness)/(i in the yxF direction 
SHPARg = (lower brick thickness)/(i in the ztf direction 
SHPAR7 = periodicity/d in the yxF direction 
SHPARg = periodicity/d in the ztf direction 



The actual numbers of dipoles A^i^, Niy, Niz, along each dimension of the upper brick, and N2X, 
N2y, N2z along each dimension of the lower brick, just be integers. Usually, A^i^^ = nint(SHPARi), 
Niy = nint(SHPAR2), A^i^ = nint(SHPAR3), N2X = nint(SHPAR4), N2y = nint(SHPAR5), N2Z = 
nint(SHPAR6), where nint(a;) is the integer nearest to x, but under some circumstances Ni^, Niy, Niz, 
N2x^ N2y, ^2z might be larger or smaller by 1 unit. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameters SHPARi 
- SHPARg: 

N = (iVi, x Niy X iVi,) + [N2x X N2y X 7V2.) (39) 
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The periodicity in the yxF and ztf directions is determined by parameters SHPARy and SHPARg. 
The physical size of the TUC is specified by the value of Ooff (in physical units, e.g. cm), specified in the 
file ddscat . par: 

d = (47r/3iV)i/3aeff (40) 

The target is a periodic structure, of infinite extent in the yxF and ztf directions. The target axes are 
set to ai = xtf - i e., normal to the "slab" - and 3.2 = yxF- The orientation of the incident radiation 
relative to the target is, as for all other targets, set by the usual orientation angles /3, 8, and $ (see fJTT] 
above) specifying the orientation of the target axes ai and a2 relative to the direction of incidence; for 
example, 8 = would be for radiation incident normal to the slab. 

The scattering directions are specified by specifying the diffraction order {M, N); for each diffrac- 
tion order one transmitted wave direction and one reflected wave direction will be calculated, with the 
dimensionless 4x4 scattering matrix S^^'^'' calculated for each scattering direction. At large distances 
from the infinite slab, the scattered Stokes vector in the (M,N) diffraction order is 

4 

IscaAM,N)=Y.^^^^^"^^^ (41) 
J = l 

where lin.j is the incident Stokes vector See Draine & Flatau (2008) for interpretation of the Sij as 
transmission and reflection efficiencies. 

20.9 SPHRN_PBC = Target consisting of group of N spheres, extended in target 
y and z directions using periodic boundary conditions 

This option causes DDSCAT to create a target consisting of a periodic array of iV— sphere structures, 
where one A^— sphere structure consists of A^ spheres, just as for target option NANSPH (see ^119. 22b . 
Each sphere can be of arbitary composition, and can be anisotropic if desired. Information for the de- 
scription of one A^-sphere structure is supplied via an external file, just as for target option NANSPH - 
see ^19^ . 

Let us refer to a single A^-sphere structure as the Target Unit Cell (TUC). The TUC is then repeated 
in the utf- and zxF-directions, with periodicities PYAEFFxceff and PZAEFFxocff, where Ocff = 
(3VTUc/47r)^/'^, where Vruc is the total volume of solid material in one TUC. This option requires 
3 shape parameters: 

DIAMX = maximum extent of target in the target frame x direction/d 
PYAEFF = periodicity in target y direction/oeff 
PZAEFF = periodicity in target z direction/aeff. 
The pertinent line in ddscat . par should read 

SHPARi SHPAR2 S HP AR3 '^^/ename' (quotes must be used) 

where 

SHPARi = target diameter in x direction (in Target Frame) in units of d 

SHPAR2= PYAEFF 
SHPAR3= PZAEFF. 

filename is the name of the file specifying the locations and relative sizes of the spheres. 

The overall size of the TUC (in terms of numbers of dipoles) is determined by parameter SHPARi, 
which is the extent of the multisphere target in the a;-direction, in units of the lattice spacing d. The 
physical size of the TUC is specified by the value of Ocff (in physical units, e.g. cm), specified in the file 
ddscat . par. 

The location of the spheres in the TUC, and their composition, is specified in file filename'. Please 
consult i il9.22l above for detailed information concerning the information in this file, and its arrangement. 

Note that while the spheres can be anisotropic and of differing composition, they can of course also 
be isotropic and of a single composition, in which case the relevant lines in the file filename ' should 
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read 

yi zi ri 1 1 1 
a;2 y2 ^2 r-2 1 1 1 
X3 ^3 rg 1 1 1 

i.e., every sphere has isotropic composition IC0MP = 1 and the three dielectric function orientation angles 
are set to zero. 

When the user uses target option SPHRN_PBC, the target now becomes a periodic structure, of infinite 
extent in the target y- and z- directions. The target axis ai ~ xlf = (I7 0, 0)tf - i-e., normal to the 
"slab" - and target axis 3.2 = ylf = (0, 1, 0)tf- The orientation of the incident radiation relative to the 
target is, as for all other targets, set by the usual orientation angles f3, Q, and <I> (see fJTT] above). The 
scattering directions are, just as for other targets, determined by the scattering angles 9s, (ps (see [j2T| 
below). 

The scattering problem for this infinite structure, assumed to be illuminated by an incident monochro- 
matic plane wave, is essentially solved "exactly", in the sense that the electric polarization of each of the 
constituent dipoles is due to the electric field produced by the incident plane wave plus all of the other 
dipoles in the infinite target. 

However, the assumed target will, of course, act as a perfect diffraction grating if the scattered radia- 
tion is calculated as the coherent sum of all the oscillating dipoles in this periodic structure: the far-field 
scattered intensity would be zero in all directions except those where the Bragg scattering condition is 
satisfied, and in those directions the far-field scattering intensity would be infinite. 

To suppress this singular behavior, we calculate the far-field scattered intensity as though the sepa- 
rate TUCs scatter incoherently. The scattering efficiency Qsca. and the absorption efficiency Qabs are 
defined to be the scattering and absorption cross section per TUC, divided by ira^g, where aoff = 
(3Vtuc /47r)^/'^, where Vtuc is the volume of solid material per TUC. 

Note: the user is allowed to set the target periodicity in the target ?/tf (or ztf) direction to values that 
could be smaller than the total extent of one TUC in the target i/tf (or ztf) direction. This is physically 
allowable, provided that the spheres from one TUC do not overlap with the spheres from neighboring 
TUCs. Note that DDSCAT does not check for such overlap. 

20.10 TRILYRPBC = Three stacked rectangular blocks, repeated periodically 

The target unit cell (TUC) consists of a stack of 3 rectangular blocks with centers on the xtf axis. The 
TUC is repeated in either the yxF direction, the ztf direction, or both. A total of 11 shape parameters 
must be specified: 

SHPARi = x-thickness of upper layer/d [material 1] 

SHPAR2 = y-width/d of upper layer 

SHPAR3 = z-width/ii of upper layer 

SHPAR4 = x-thickness/d of middle layer [material 2] 

SHPAR5 = y-width/d of middle layer 

SHPARg = z-width/rf of middle layer 

SHPAR7 = x-width/d of lower layer 

SHPARg = y-width/d of lower layer 

SHPARg = z-width/rf of lower layer 

SHPARio = period/rf in y direction 

SHPARii= period/d in z direction 

21 Scattering Directions 
21.1 Isolated Finite Targets 

DDSCAT calculates scattering in selected directions, and elements of the scattering matrix are reported 
in the output files vxxxryy'Kzzz . sea . The scattering direction is specified through angles 9s and 0s (not 
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to be confused with the angles Q and $ which specify the orientation of the target relative to the incident 
radiation!). 

For isolated finite targets [i.e., PBC not employed) there are two options for specifying the scatter- 
ing direction, with the option determined by the value of the string CMDFRM read from the input file 

ddscat . par. 

1. If the user specifies CMDFRM=' LFRAME' , then the angles 9, (j) input from ddscat .par are 
understood to specify the scattering directions relative to the Lab Frame (the frame where the 
incident beam is in the x— direction). 

When CMDFRM=' LFRAME', the scattering angle 0s is simply the angle between the incident 
beam (along direction x) and the scattered beam {Os ~ for forward scattering, 9s — 180° for 
backscattering). 

The scattering angle 0s specifies the orientation of the "scattering plane" relative to the x, y plane 
in the Lab Frame. When 0^ = the scattering plane is assumed to coincide with the x, y plane. 
When 0s = 90° the scattering plane is assumed to coincide with the x, z plane. Within the 
scattering plane the scattering directions are specified by < 9s < 180°. 

2. If the user specifies CMDFRM=' TFRAME', then the angles 9, input from ddscat .par are 
understood to specify the scattering directions fig relative to the Target Frame (the frame defined 
by target axes ai, a2, §.3). 8 is the angle between ai and fig, and is the angle between the ai — n 
plane and the ai — 3.2 plane. 

Scattering directions for which the scattering properties are to be calculated are set in the parameter 
file ddscat . par by specifying one or more scattering planes (determined by the value of 0s) and for 
each scattering plane, the number and range of 9s values. The only limitation is that the number of 
scattering directions not exceed the parameter MXSCA in DDSCAT . f (in the code as distributed it is set 

toMXSCA=1000). 



21.2 Scattering Directions for Targets that are Periodic in 1 Dimension 

For targets that are periodic, scattering is only allowed in certain directions (see Draine &. Flatau 2008). 
If the user has chosen a PBC target (e.g, NSPPBC, CYLPBC, HEXPBC, or RCTPBC) but has set one of 
the periodicities to zero, then the target is periodic in only one dimension - e.g., CYLPBC could be used 
to construct an infinite cylinder 

In this case, the scattering directions are specified by giving an integral diffraction order M = 
0, ±1, ±2, ... and one angle, the azimuthal angle around the target symmetry axis. The diffraction 
order M determines the projection of ks onto the "repeat" direction. 

For example, if PYD > (target repeating in the ?;tf direction), then M determines the value of 
ksy = kg • yxF, where yxp is the unit vector in the Target Frame y— direction; 

ksy = koy + 2TTM/Ly (42) 

where koy = ko • yxF, where ko is the incident k vector. Note that the diffraction order M must satisfy 
the condition 

{koy - ko){Ly/2Tr) < M < (fco - koy){Ly/2Tr) (43) 

where Ly = PYD x d is the periodicity along the y axis in the Target Frame. i\/ = is always an allowed 
diffraction order. 

The azimuthal angle ^ defines a right-handed rotation of the scattering direction around the target 
repetition axis. Thus for a target repetition axis yxF, 

ksx = /c^cosC , (44) 
ksz = k±sm( , (45) 

where k± = (fcg — k'^y)^^'^, with ksy = koy + 2'kM/ Ly. For a target with repetition axis ztf, 

ksx = fc^cosC , (46) 
ksy = /c^sinC , (47) 
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where fcj_ = {B - k1,f'^, k,, = ko, + 27:M/L,. 

The user selects a diffraction order M and the azimuthal angles ( to be used for that M via one 
line in ddscat . par. An example would be to use CYLPBC to construct an infinite cylinder with the 
cylinder direction in the yxF direction: 

' =================== Parameter file ===================' 

'**** PRELIMINARIES ****' 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 
' **** Initial Memory Allocation 

50 50 50 = size allowance for target generation 
' **** Target Geometry and Composition 
'CYLNDRPBC' = CSHAPE*9 shape directive 

2.0 20. 2 2.0 0. = shape parameters SHPARl, SHPAR2, SHPAR3, SHPAR4, SHPAR5, 

1 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/ml . 33_0 . 01' = name of file containing dielectric function 
' -k-k-k-k Error Tolerance 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC| E>-ACA | X> ) / (NORM OF AC|E>) 

'-k-k-k-k Interaction cutoff parameter for PBC calculations -kk-k-k' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'-k-k-k-k Angular resolution for calculation of <cos>, etc. 

0.5 = ETASCA (number of angles is proportional to [ (3+x) /ETASCA] "2 ) 

'-k-k-k-k Wavelengths (micron) 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
' -k-k-k-k Effective Radii (micron) -k-k-k-k ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
'-k-k-k-k Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient.) 

I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) ) 
' kkkk Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

60. 60. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, TRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first TORI (0 to begin fresh) 
'-k-k-k-k Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 
' **** Specify Scattered Directions ****' 

'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 

0. 0. 350. 30 = M, zeta_min, zeta_max, dzeta (in degrees) 

In this example, a single diffraction order M = is selected, and C is to run from C,-a\\\\ = to 
Cmax = 360° in increments of = 30°. 

There may be additional lines, one per diffraction order. Remember, however, that every diffraction 
order must satisfy eq. (03]). 

21.3 Scattering Directions for Targets for Doubly-Periodic Targets 

If the user has specified nonzero periodicity in both the y and z directions, then the scattering directions 
are specified by two integers - the diffraction orders M, for the y, z directions. The integers M and 
iV must together satisfy the inequality 

(k^y + 2^M/Lyf + (fco. + 2^N/L,f < kl (48) 

which, for small values of Ly and L^, may limit the scattering to only (il/, N) ~ (0,0) [(0,0) is always 
allowed). 

DDSCAT will automatically calculate forward scattering for each selected {M,N); this will be 

followed by calculation of backscattering for the same set of (A/, A^) values. 
Here is a sample ddscat . par file as an example: 



INCIDENT POLARIZATION STATE 



' =================== Parameter file ' 

PRELIMINARIES ****' 

'NOTORQ' = CMT0RQ*6 (DOTORQ, NOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) 
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR) 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 
' -k-k-k-k Initial Memory Allocation 

50 50 50 = dimensioning allowance for target generation 
'k-k-k-k Target Geometry and Composition 
'RCTGL_PBC' = CSHAPE*9 shape dirctive 

20.0 1.0 1.0 1.0 1.0 = shape parameters PARI, PAR2, PAR3,... 

1 = NCOMP = number of dielectric materials 

' /u/draine/work/DDA/diel/ml . 33_0 . 01' = name of file containing dielectric function 
CONJUGATE GRADIENT DEFINITIONS 

= INIT (TO BEGIN WITH | X0> = 0) 

l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC I E>-ACA | X> ) / (NORM OF AC|E>) 

' kkkk Interaction cutoff parameter for PBC calculations kkkk' 

5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 

'kkkk Angular resolution for calculation of <cos>, etc. -k-k-k-k' 

0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] " 2 ) 

'kkkk Wavelengths (micron) 

6.283185 6.283185 1 ' INV = wavelengths (first, last, how many, how=LIN, INV, LOG) 
' k-k-k-k Effective Radii (micron) -k-k-k-k ' 

2 2 1 'LIN' = eff. radii (first, last, how many, how=LIN, INV, LOG) 
' **** Define Incident Polarizations 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient.) 

I = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) ) 
'-k-k-k-k Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

60. 60. 1 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 
'kkkk Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 

II 12 21 22 31 41 = indices ij of elements to print 

Specify Scattered Directions ****' 
'TFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

1 = number of scattering orders 
= M, N (diffraction orders) 



22 Incident Polarization State 

Recall that the "Lab Frame" is defined such that the incident radiation is propagating along the xlf 
axis. DDSCAT allows the user to specify a general elliptical polarization for the incident radiation, by 
specifying the (complex) polarization vector Gqi. The orthonormal polarization state eo2 = xlf x 
is generated automatically if ddscat .par specifies I0RTH=2. 

For incident linear polarization, one can simply set Gqi = y by specifying 
(0,0) (1,0) (0,0) 

in ddscat . par; then eo2 = z- For polarization mode Gqi to correspond to right-handed circular po- 
larization, set eoi = (y + by specifying (0,0) (1,0) (0,1) in ddscat . par (DDSCAT 
automatically takes care of the normalization of Gqi); then eo2 = («y + z)/^2, corresponding to left- 
handed circular polarization. 
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56 




30 60 90 120 150 180 



Figure 8: Scattered intensity for incident unpolarized light for a cube with m = 1.5 + O.OOli and x — 27raoft/A = 6.5 (i.e., 
d/X — 1.6676, where d is the length of a side). The cube is tilted with respect to the incident radiation, with O = 30°, and rotated 
by /3 = 15° around its axis to break reflection symmetry. The Mueller matrix element Sn is shown for 4 different scattering planes. 
The strong forward scattering lobe is evident. It is also seen that the scattered intensity is a strong function of scattering angle (ps as 
well as 9s - at a given value of 63 (e.g., 9s = 120°), the scattered intensity can vary by orders of magnitude as (j>s changes by 90°. 



23 Averaging over Scattering Directions: ^'(1) = (cos 6*5), etc. 
23.1 Angular Averaging 

An example of scattering by a nonspherical target is shown in Fig. [H showing the scattering for a tilted 
cube. Results are shown for four different scattering planes. 

DDSCAT automatically carries out numerical integration of various scattering properties, including 

• (cos6's); 

• g = (cos6's)xlf + (sin 0s cos(j>s)y + (sin 0s siii0s)z (see flSjl; 

• Qr, provided option DOTORQ is specified (see i fTsT l. 

The angular averages are accomplished by evaluating the scattered intensity for selected scattering di- 
rections {ds, 4>s), and taking the appropriately weighted sum. Suppose that we have Ng different values 

of 0s, 

= 0,, j = l,...,Ng , (49) 
and for each value of 0j, N^{j) different values of (j)s: 

(/)s=0,,fc, k = l,...,N4j) . (50) 

For a given j, the values of (j)j^k are assumed to be uniformly spaced: (j)j,k+i — <t>j,k = '2TT/N^{j). The 
angular average of a quantity f{0s, 4>s) is approximated by 

{f)^Y sin0sd0s / d0s/(0s,0s) ~ ^Y^Y. (51) 
= 77T^[cos(^J-i)~cos(0,+i)] , j=2,...,Ng-l . (52) 
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^l,k 



2-K 
27r 



COs(0i) + COS(02) 



C0s(6'Ar,_i) +C0s(6'ArJ 



(53) 
(54) 



N^{Ng) 

For a sufficiently large number of scattering directions 

Ne 

Nsc. = Y,N4j) (55) 

the sum ( ISTT i approaches the desired integral, but the calculations can be a significant cpu-time burden, 
so efficiency is an important consideration. 

23.2 Selection of Scattering Angles Og , (ps 

Beginning with DDSCAT 6.1, an improved approach is taken to evaluation of the angular average. Since 
targets with large values of x = 27racff / A in general have strong forward scattering, it is important to 
obtain good sampling of the forward scattering direction. To implement a preferential sampling of the 
forward scattering directions, the scattering directions {6, (f>) are chosen so that 9 values correspond 
to equal intervals in a monotonically increasing function s{9). The function s{6) is chosen to have a 
negative second derivative cPs/dO'^ < so that the density of scattering directions will be higher for 
small values of 9. DDSCAT takes 

■K^) = ^ + 7r^ • (56) 



This provides increased resolution (i.e., increased ds/dO) for 9 < 
scattering lobe, so we take 



We want this for the forward 



The 9 values run from 6' = to 6* 



1 + x 

TT, corresponding to uniform increments 
1 



As 



Ng-1 



If we now require that 

max[A0 

we determine the number of values of 9: 

Ne = l + 



As 



tt/2 



{ds/d9)g=Tr '3 + x 
2(3 + a;) [1 + l/(7r 



(57) 



(58) 



(59) 



(60) 



77 [l + 0o/(^ + 0o)2] ■ 

Thus for small values of x, max[A6'] = 30°7], and for a; ^ 1, max[A0] 90°?7/a;. For a sphere, 
minima in the scattering pattern are separated by ^ 180°/a:;, so 77 = 1 would be expected to marginally 
resolve structure in the scattering function. Smaller values of rj will obviously lead to improved sampling 
of the scattering function, and more accurate angular averages. 

The scattering angles 9j used for the angular averaging are then given by 

9o) + [(1 + 9o- sjf + i9os,] 



where 



U - 1)^ 



For each ( 



Ng-1 

we must choose values of ( 



1 



2 
1 



J = h...,Ng 



(61) 



(62) 



For 9i = Q and 9^^ = tt only a single value of cj) is needed 



(the scattering is independent of in these two directions). For < 9j < tt we use 

= max{3,nint [ii: sm{9j)/ {9j+i - 9j^i)]} , 
where nint = nearest integer. This provides sampling in consistent with the sampling in 9. 



(63) 
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Figure 9: Errors in (cos 6) and (cos'^ 9) calculated for z N — 17904 dipole pseudosphere with m — 1.33 + O.Oli, as functions 
of a:: = 2na/X. Results are shown for different values of the parameter 77. 
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Figure 10: Same as Fig.[9l but f or a Af = 32768 dipole cube with m = 1.5 + O.OOli. 



24 SCATTERING B Y FINITE TARGETS: THE MUELLER MATRIX 



23.3 Accuracy of Angular Averaging as a Function of 77 

Figure |9] shows the absolute errors in {cos 9) and (cos^ 9) calculated for a sphere with refractive index 
771 = 1.33 + O.Oli using the above prescription for choosing scattering angles. The error is shown as a 
function of scattering parameter x. We see that accuracies of order 0.01 are attained with 1] = 1, and that 
the above prescription provides an accuracy which is approximately independent of x. We recommend 
using values of 77 < 1 unless accuracy in the angular averages is not important. 



24 Scattering by Finite Targets: The Mueller Matrix 
24.1 Two Orthogonal Incident Polarizations (I0RTH=2) 

Throughout the following discussion, ic, y, z are unit vector defining the Lab Frame (thus ko = x). 

DDSCAT internally computes the scattering properties of the dipole array in terms of a complex 
scattering matrix fmi{ds, 4's) (Draine 1988), where index 1 — 1,2 denotes the incident polarization state, 
m = 1, 2 denotes the scattered polarization state, and 9s,(j>s specify the scattering direction. Normally 
DDSCAT is used with I0RTH=2 in ddscat .par, so that the scattering problem will be solved for 
both incident polarization states (1 = 1 and 2); in this subsection it will be assumed that this is the case. 

Incident polarization states ^ = 1, 2 corresp ond t o polarization states eoi, eo2; recall that polarization 
state eoi is user-specified, and eo2 = x x Sqi '1'"* Scattered polarization state m = 1 corresponds to 



linear polarization of the scattered wave parallel to the scattering plane (ei = ep^ = 9g) and m = 2 

corresponds to linear polarization perpendicular to the scattering plane (in the +<j>s direction). The 
scattering matrix /,„/ was defined (Draine 1988) so that the scattered electric field Es is related to the 
incident electric field (0) at the origin (where the target is assumed to be located) by 



E, • 9s \ ^ cxp(ik, ■ r) / /n f,^ \ f E,(0) • e*, 
E, -0, J kr \ /21 /22 / I E,(0)-e52 



(64) 



The 2x2 complex scattering amplitude matrix (with elements 5*1, 6*2, 5*3, and 5*4) is defined so that (see 
Bohren & Huffman 1983) 



E, • \ cxp(ik, ■ r) f S2 S3\ f E,(0) • e. 



E,, • 0,, / -ikr V -^4 Si J \ E,(0) • e,j_ 



(65) 



where e^n , e^j^ are (real) unit vectors for incident polarization parallel and perpendicular to the scattering 
plane (with the customary definition of e^^ = 6^11 x xlf)- 



From ( 164165b we may write 



^2 ^3 W E,(0).e,|| A _Y '^12 A f E,(0).e5i 

Si Si ){ E,(0)-e,i J -hi -/22 ) \ E.(0)-e52 



(66) 



Let 



(67) 
(68) 
(69) 
(70) 

Note that since eoi, eo2 could be complex (i.e., elUptical polarization), the quantities a, 6, c, d are com- 
plex. Then 



a — 


^01 


YLF 


b = 


^01 


zlf 


c = 


p* 
'='02 


yLF 


d = 


^02 


zlf 



Draine (1988) adopted the convention eo2 = x x eSi. The customary definition of ei± is e^^ = ei|| x k. Thus, if eoi 
then eo2 = —ei±. 

A frequent choice is eoi = Ylf, with eo2 = xlf x ylf ~ zlf- 
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and eq. i66i can be written 

^2 ^3 \ / E,(0) -6,11 \ / -/n -/i2 \f a b\f E,(0) -YLF 

^4 Si ){ E,(0)-e,i ) \ f2l f22 )\C d ) \ E,(0)-ZLF 

The incident polarization states ej|| and e^j^ are related to xlf, Ylf, zlf by 

Ylf \ ^ ( cos<j>s sin \ / e^n 
ZLF / \ sin0^. -cos 0s J \ ei± 



(72) 



(73) 



substituting ( l73b into ( |72] | we obtain 

^2 ^3 W E,(0) • e,|| \ / -/ii -/i2 5 \ / COS0,, sin0,, \ / E,(0) • e,,, 

^4 Si )\ E,(0) • y * 1^ /21 f22 ) \ c d){ sin0, -cos0, E,(0) • e,^ 

(74) 

Eq. ( |74l i must be true for all Ei(0), so we obtain an expression for the complex scattering amplitude 
matrix in terms of the : 



/ ^2 


S3 ^ 




I ^4 


Si ) 


-A 



-/ii -/12 \(^ & W cos 0s sin0, 
/21 f22 J \ c d J \ sin 0s -cos( 



(75) 



This provides the 4 equations used in subroutine GETMUELLER to compute the scattering amplitude 
matrix elements: 

51 = -i [/2i(6cos0s - asin0s) + /22(rfcos0s - csin0s)] , (76) 

52 = [/ii(acos0s + fesin0s) + /i2(ccos0s + rfsin0s)] , (77) 

53 = i [/ii(5cos0s - asin0s) + /i2(rfcos0s - csin0s)] , (78) 
<S'4 = i [/2i(acos0s + ^^sin0s) + /22(ccos0s + cfsin0s)] • (79) 



24.2 Stokes Parameters 

It is both convenient and customary to characterize both incident and scattered radiation by 4 "Stokes 
parameters" - the elements of the "Stokes vector". There are different conventions in the literature; we 
adhere to the definitions of the Stokes vector (I,Q,U,V) adopted in the excellent treatise by Bohren & 
Huffman (1983), to which the reader is referred for further detail. Here are some examples of Stokes 
vectors (/, Q, [/, V) = (1, Q//, U/I, V/I)I: 

• (1,0, 0,0)/: unpolarized light (with intensity /); 

• (1,1, 0, 0)/ : 100% linearly polarized with E parallel to the scattering plane; 

• (1,-1,0, 0)/ : 100% Unearly polarized with E perpendicular to the scattering plane; 

• (1, 0, 1, 0)/ : 100% linearly polarized with E at +45° relative to the scattering plane; 

• (1,0,-1, 0)/ : 100% linearly polarized with E at -45° relative to the scattering plane; 

• (1, 0, 0, 1)/ : 100% right circular polarization {i.e., negative helicity); 

• (1,0, 0, —1)/ : 100% left circular polarization (i.e., positive helicity). 



24.3 Relation Between Stokes Parameters of Incident and Scattered Radiation: 
The Mueller Matrix 



It is convenient to describe the scattering properties in terms of the 4x4 Mueller matrix relating the 





K)and (/s, 


/ /. \ 


/ 


Qs 


1 


Us 


^.2y,2 


V Vs J 


V 



Sii Si2 Si3 Si4 
S2I S22 S23 S24 
S31 S32 S33 S34 



I \ 



(80) 
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Once the amplitude scattering matrix elements are obtained, the Mueller matrix elements can be com- 
puted (Bohren & Huffman 1983): 





= (|^1P + |^2P + |^3P + 


1^4^) 


/2 


Sl2 


= (|52p-|5lP + |54p- 


1^3^ 


/2 


5*13 


= RC(52^3*+^1^4*) , 






5*14 


= ImiS2S*3- SiSl) , 






5*21 


= {\S2\^^\Si\' + \S3\'- 




/2 


S22 


= {\Sl\^ + \S2\'~\S3\^^ 


1^4^ 


/2 


S23 


= Rc{S2S;~ SiSl) , 






S24 


= lm{S2S*3+SiSl) , 






S31 


= RciS2Sl + SiS;) , 






S32 


= Rc{S2Sl~SiS;) , 






S33 


= ReiSiS^ + SsSl) , 






S34 


= Im{S2S*,+ SiS;) , 






S41 


= Im (545* + 5*15*) , 






S42 


= ImiSiS;- SiS^) , 






843 


= lm{SiS*2- S3SI) , 






5*44 


= Rc (5*1 5*2 — 53^4) 







These matrix elements are computed in DDSCAT and passed to subroutine WRITESCA which handles 
output of scattering properties. Although the Muller matrix has 16 elements, only 9 are independent. 

The user can select up to 9 distinct Muller matrix elements to be printed out in the output files 
wxxxryykzzz . sea and wxxxryyori . avg; this choice is made by providing a list of indices in ddscat . par 
(see Appendix IaIi. 

If the user does not provide a list of elements, WRITESCA will provide a "default" set of 6 selected 
elements: ^21, 5*31, 6*41 (these 4 elements describe the intensity and polarization state for scattering 
of unpolarized incident radiation), 5i2, and £'13. 

In addition, WRITESCA writes out the linear polarization P of the scattered light for incident unpo- 
larized light (see Bohren & Huffman 1983): 



24.4 Polarization Properties of the Scattered Radiation 



The scattered radiation is fully characterized by its Stokes vector {Is,Qs,Us,Vs). As discussed in 
Bohren & Huffman (eq. 2.87), one can determine the linear polarization of the Stokes vector by op- 
erating on it by the Mueller matrix of an ideal linear polarizer: 



/ 1 



Spol 



cos 2^ 
cos 2^ cos^ 2^ 
sin 2^ sin 2^ cos 2^ 




sin 2^ 
cos 2^ sin 2^ 
sin^ 2^ 






/ 



(83) 



where ^ is the angle between the unit vector 6s parallel to the scattering plane ("SP") and the "transmis- 
sion" axis of the linear polarizer Therefore the intensity of light polarized parallel to the scattering plane 
is obtained by taking f = and operating on {Is,Qs,Us,Vs) to obtain 



SP) = lih 



Qs) 



(84) 



2fc2r^ 



+ S2l)h + (5i2 + S22)Q^ + {Sl3 + S23)U, + (5i4 + ^24)^,] . (85) 
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Similarly, the intensity of light polarized perpendicular to the scattering plane is obtained by taking 

? = ^2: 

/(E, ±SP) = lils-Qs) (86) 

= - S2l)h + iSl2 - S22)Q^ + [Su - S2z)U, + [Su " &4)F,] . (87) 



2fc2r 

24.5 Relation Between Mueller Matrix and Scattering Cross Sections 

Differential scattering cross sections can be obtained directly from the Mueller matrix elements by noting 
that 

l_ f dCs, 



Is = -[ ] h ■ (88) 



Let SP be the "scattering plane": the plane containing the incident and scattered directions of propaga- 
tion. Here we consider some special cases: 

• Incident light unpolarized: Stokes vector Si = /(I, 0, 0, 0): 

- cross section for scattering with polarization Es || SP: 

^ = ^ (l^2p + l^sP) ^ ^ (5n + S2.) (89) 

- cross section for scattering with polarization E^ _L SP: 

^ = ^(l^iP + l^4n = ^(^n-^2i) (90) 

- total intensity of scattered light: 

^ = i (l^iP + \S2? + I^3P + \S,n = ^^n (91) 

• Incident light polarized with E^ || SP: Stokes vector Si ~ /(1, 1, 0, 0): 

- cross section for scattering with polarization Eg || to SP: 



cross section for scattering with polarization E^ -L to SP 



7^1^21' = :.^ (511+^12+521 +522) (92) 



dCsca _ _1_ I o |2 _ _j_ 

d^l ^ pl*^*! - 2fc2 



T7l^4|' = x77(5ii + 512-521-522) (93) 



total scattering cross section: 



^(|52p + |54n = ^(5ii + 5i2) (94) 



dVL F ^' ' 1 -1 / ^, 

Incident light polarized with Ej + SP: Stokes vector Si = /(I, —1, 0, 0): 
- cross section for scattering with polarization Eg || to SP: 



- cross section for scattering with polarization E^ + SP: 
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- total scattering cross section: 



"^^^^^ ^ {\S,\' + \S3n = ^{Su-S,2) (97) 



• Incident light linearly polarized at angle 7 to SP: Stokes vector = /(I, cos 27, sin 27, 0): 
- cross section for scattering with polarization Eg || to SP; 

dCsca, 1 



7777 [(^11 + ^21) + (^12 + S22) COS 27 + (5i3 + ^23) sin 27] (98) 



dn 2fc2 

- cross section for scattering with polarization -L SP: 

dCsca, 1 



dn 2fc2 

- total scattering cross section: 

dCsca 1 



7772 [(^11 - ^21) + (^12 - ^22) COS 27 + (5i3 - 523) sill 27] (99) 



dn fc2 



[S-ii +S'i2Cos27 + 5i3sin27] (100) 



24.6 One Incident Polarization State Only (I ORTH= 1) 

In some cases it may be desirable to limit the calculations to a single incident polarization state - for 
example, when each solution is very time-consuming, and the target is known to have some symmetry 
so that solving for a single incident polarization state may be sufficient for the required purpose. In this 
case, set I0RTH=1 in ddscat .par. 

When I0RTH=1, only /n and /21 are available; hence, DDSCAT cannot automatically generate 
the Mueller matrix elements. In this case, the output routine WRITE SCA writes out the quantities |/ii p, 
1/21 Re(/ii/|i), and Im(/ii/|j) for each of the scattering directions. 

The differential scattering cross section for scattering with polarization Es || and _L to the scattering 
plane are 

/ ^r' \ 1 

l/iiP (101) 
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V dn 


fc2 


( dCsca. \ 


1 




fc2 



f2i|' (102) 



Note, however, that if IPHI is greater than 1, DDSCAT will automatically set I0RTH=2 even if 
ddscat .par specified I0RTH=1: this is because when more than one value of the target orientation 
angle $ is required, there is no additional "cost" to solve the scattering problem for the second incident 
polarization state, since when solutions are available for two orthogonal states for some particular target 
orientation, the solution may be obtained for another target orientation differing only in the value of 
$ by appropriate linear combinations of these solutions. Hence we may as well solve the "complete" 
scattering problem so that we can compute the complete Mueller matrix. 



25 Scattering by Periodic Targets: Generalized Mueller Matrix 

The Mueller scattering matrix formalism (e.g., Bohren & Huffman 1983) was originally defined to de- 
scribe scattering of incident plane waves by finite targets, such as aerosol particles or dust grains. 

Draine & Flatau (2008) have extended the Mueller scattering matrix formalism to also apply to 
targets that are periodic and infinite in one or two dimensions (e.g., an infinite chain of particles, or a 
two-dimensional array of particles). 

The scattering matrix reported by DDSCAT in the vaaarbbb'kccc . sea output files corre- 

spond to the ordinary Mueller scattering matrix for single finite targets (periodic in n = dimensions), 
or to the scattering matrices S^^'^'^ for targets periodic in one dimension, or to the scattering matrices 
S^'^'^^ for scattering by two-dimensional arrays. See Draine & Flatau (2008) for the exact definitions of 



COMPOSITE TARGETS WITH ANISOTROPIC CONSTITUENTS 



26 Composite Targets with Anisotropic Constituents 

Section[T9]includes targets composed of anisotropic materials (ANIELLIPS, ANIRCTNGL, ANI_ELLJ2, 
ANI_ELL_3). However, in each of these cases it is assumed that the dielectric tensor of the target ma- 
terial is diagonal in the "Target Frame". For targets consisting of a single material (in a single domain), 
it is obviously possible to choose the "Target Frame" to coincide with a frame in which the dielectric 
tensor is diagonal. 

However, for inhomogeneous targets containing anisotropic materials as, e.g., inclusions, the optical 
axes of the constituent material may be oriented differently in different parts of the target. In this case, it 
is obviously not possible to choose a single reference frame such that the dielectric tensor is diagonalized 
for all of the target material. 

To extend DDSCAT to cover this case, we must allow for off-diagonal elements of the dielectric 
tensor, and therefore of the dipole polarizabilities. It will be assumed that the dielectric tensor e is 
symmetric (this excludes magnetooptical materials - check this). 

Let the material at a given location in the grain have a dielectric tensor with complex eigenvalues e*-^-* , 
j = 1 — 3. These are the diagonal elements of the dielectric tensor in a frame where it is diagonalized 
(i.e., the frame coinciding with the principal axes of the dielectric tensor). Let this frame in which the 
dielectric tensor is diagonalized have unit vectors corresponding to the principal axes of the dielectric 
tensor. We need to describe the orientation of the "Dielectric Frame" (DF) - defined by the "dielectric 
axes" ej - relative to the Target Frame. It is convenient to do so with rotation angles analogous to the 
rotation angles O, (f>, and (3 used to describe the orientation of the Target Frame in the Lab Frame: 
Suppose that we start with unit vectors aligned with the target frame Xj . 

1. Rotate the DF through an angle ^df around axis yxF, so that Ouf is now the angle between ei 
and XTF . 

2. Now rotate the DF through an angle (p^p around axis xtfj in such a way that 62 remains in the 
Xtf ^ Gi plane. 

3. Finally, rotate the DF through an angle /3df around axis ei. 

The unit vectors are related to the TF basis vectors Xtf, Ytf- ^tf by: 

ei = Xtf cos 6'df + Ytf sin 6'df cos 0df + ztf sin 6'df sin (/)df (103) 

62 = —Xtf sin 6'df cos /3df + Ytf [cos 6'df cos /3df cos (/)df — sin /3df sin 0df] 

+ztf[cos0df cos/3df sin(/)DF + sin/3DF cos^df] (104) 

63 = Xtf sin 6'df sin /3df — Ytf [cos ^df sin /3df cos 0df + cos /3df sin (/)df] 

— ztf[cos 9y)f sin /3df sin 0df — cos /3df cos c/)df] (105) 

or, equivalently: 

Xtf = ei cos 6'df — ^2 sin 6'df cos /?df + e3 sin 6'df sin /3df (106) 
Ytf = Gi sin 6'df cos 4>df + ^2 [cos 9uf cos /3df cos (pBF — sin /?df sin 4>f)f] 

— e3[cos0DF sin/9DF cos^DF + cos /3df sin (/)df] (107) 
ztf = ei sin ^df sin (/)df + ^2 [cos Obf cos /3df sin ^df + sin /3df cos 0df] 

—63 [cos 6'df sin /3df sin (puF — cos Pbf cos (pup] (108) 

Define the rotation matrix Rij = x^ • : 

R^J = (109) 

fcos duF sin Oof cos (j!>df sin ^df sin 0df \ 

-sin ^F cos /3df cos Sdf cos /9df cos 0df — sin /3df sin <j!>df cos Sdf cos /9df sin ^DF+sin /3df cos (j!>df 1 
sin 9uF sin /3df —cos ^f sin /3df cos ^df — cos Puf sin 0df —cos ^df sin /3df sin (/)df + cos /3df cos <I)df/ 

(110) 

and its inverse 

{R'')ij= (111) 
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(cos 9uF — sin 6'df cos /?df sin 9df sin /3df \ 

sin duF cos 0DF cos 6'df cos Puf cos 0DF ~ sin /3df sin 0df — cos 6'df sin /3df cos c^df — cos /3df sin (j!)df 1 
sin doF sin (;/>df cos 6uf cos /3df sin ^df + sin /3df cos ^df — cos Oof sin /3df sin ^df + cos /3df cos 4>f>fJ 

(112) 

The dielectric tensor is diagonal in the DF. If we calculate the dipole polarizabihlity tensor in the DF it 
will also be diagonal, with elements 

\ 

cP^=\ a^F (113) 
The polarizability tensor in the TF is given by 

{a^'')„n = i?,,(«°^),fe(i?-')fcm (114) 



diagonal elements of the dielectric tensor in the DF, and the three rotation angles 6'df, /?df, and 
0DF- We first use the LDR prescription and the elements to obtain the three diagonal elements a^J^ of 
the polarizability tensor in the DF. We then calculate the polarizability tensor a^^ using eq. ( II 14b . 



27 Postprocessing: Calculating E and B In or Near the Target 

For some scientific applications (e.g., surface-enhanced Raman scattering) one wishes to calculate the 
electric field just outside the target boundary. It may also be of interest to calculate the electromagnetic 
field at positions within the target volume. DDSCAT includes tools to carry out such calculations. 

When DDSCAT is run with option IWRP0L=1, it creates files vixxxryyy'Kzzz . poln for «=1 and (if 
I0RTH=2) «=2. These files store the solution for the dipole polarizations, allowing postprocessing. Of 
particular interest is calculation of the electromagnetic field at locations in or near the target. 

To facilitate this, we provide a Fortran subroutine READPOL (in the file readpol . f 90) that reads 
in data from the vixxxryyy'k.zzz . poln file. 

In addition, we provide a Fortran program DDfield.f90 that calls subroutine READPOL, and uses 
the data read in by READPOL to calculate the electric and magnetic field at locations specified by the 
user. The E and B contributions from the oscillating dipoles are calculated [including retardation effects 
- see, e.g., equations (7) and (10) from Draine & Weingartner (1996)] giving the "exact" electromagnetic 
field contributed by the oscillating dipoles. The incident plane wave field is added to the field generated 
by the oscillating dipoles to obtain the total E and B at the specified location. 

The user is interested in the continuous and finite field in a real target, as opposed to the field gener- 
ated by an array of oscillating point dipoles, which diverges as one approaches any one of the dipoles. 
To suppress the divergences, the contributions of dipoles within a distance d of the selected point is ob- 
tained by multiplying the "exact" contribution from the point dipole by a factor 0(r) that goes to zero 
more rapidly than r^, thus ensuring that the electric field contribution from an individual dipole goes to 
zero as we approach the dipole location. This is desired, because the discrete dipole approximation is 
based on the assumption that a given dipole is polarized by the electric field from the incident wave plus 
the electric field contributed by all other dipoles. When dealing with periodic (i.e., infinite) targets, the 
function (f){r) smoothly suppresses the contributions from distant dipoles, allowing the summations to be 
truncated. Please refer to Draine & Flatau (2008) for further discussion of (f){r) and examples of results 
calculated byDDfield. 

To calculate the electric and magnetic field at user-selected locations: 

• Run ddscat with WRP0L=1 to create a file with the stored polarization (for example, wOOOrOOOkOOO.poll) 

• Compile and link the DDf ield code, by entering the src directory and typing 

make ddfield 

This will create an executable named DDf ield 

• Use an editor to create a file ddfield. in providing the name of the file {in quotes!) containing 
the stored target polarization. 
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followed by the value of the interaction cutoff parameter 7 (see Draine & Flatau 2008) to be used, 

followed by the locations where the E and B field is to be calculated, e.g: 

filename 

7 

xi yi zi 

X2 2/2 Z2 

Xn UN Zn 



Here Xj yj Zj are integer or decimal numbers giving the coordinates (in the "Target Frame", in 
units of the lattice spacing d) of the positions at which the electric field is to be evaluated. For 
example: 



'wOOOrOOOkOOO .poll' 
5 . OOe- 
-35 . 



name of file with stored polarization 



gamma (interaction cutoff parameter) 
. = x/d, y/d, z/d 
. 



















The above example will sample the EM field at 17 points. Note that x, y, z are specified in the 
Target Frame, and in units of d. You will need to know the location of the TF origin (0, 0, 0) - see 
the target descriptions in §|T9l20] 

• Type 

ddf ield 
to run the executable. 

• The output will be written to two files, DDf ield . E and DDBf ield . B, containing the electric 
and magnetic fields (given as complex vectors) at each chosen location. The fields are evaluated at 
a time when the factor e^"^* is 1 (e.g., f = 0, 27r/(jj, ....). 

We have run DDf ield on the output created by the "test problem" of a rectangular target with 
refractive index m = 1.33 + O.Oli (see fj8]l. 
The DDf ield. E file will look like 

12288 = number of dipoles in Target 





1 


32 = JXMIN 


JXMAX 






1 


24 = JYMIN 


JYMAX 






1 


16 = JZMIN 


JZMAX 




-33 


000000 


. 000000 


= (x_TF/d)min, {x_TF/d)max 




-12 


000000 


12.000000 


= (y_TF/d)min, (y_TF/d)max 




-8 


000000 


8.000000 


= (z_TF/d)min, (z_TF/d)max 




-4 


610482 


0.000000 


= xmin (TF) , xmax (TF) (phys. 


units) 


-1 


676539 


1 . 676539 


= ymin (TF) , ymax (TF) (phys. 


units) 


-1 


117693 


1 . 117693 


= zmin (TF) , zmax (TF) (phys. 


units ) 





139712 


= d (phys units) 







000000 


. 000000 


= PYD,PZD = period_y/dy, period_z/dz 





000000 


. 000000 


= period_y, period_2 (phys . 


units) 


6 


283185 


= wavelength in ambient medium (phys units) 





139712 


= k_x*d for 


incident wave 







000000 


= k_y*d for 


incident wave 





FINALE 



0.000000 = k_z*d for incident wave 



5 . OOOE-03 


= gamma ( 


parameter for summation cutoff) 














0.000000 





000000 


_ 


(Re, 


Im) Einc, 


X (0 


0,0) 


















1.000000 





000000 


_ 


(Re, 


Im) E_inc, 


Y (0 


0,0) 


















0.000000 





000000 




(Re, 


Im) E_inc, 


z (0 


0,0) 


















x/d 




ii/d 




z/d 




E_x 




_ 


E_ 


_y - 





_ 


E_ 


_z 




-35.0000 





0000 





. 0000 


0.00000 





00000 





18730 





79489 





00000 





00000 


-34.0000 





0000 





. 0000 


0.00000 





00000 





02038 





80327 





00000 





00000 


-33.0000 





0000 





. 0000 


0.00000 





00000 


-0 


15115 





79529 





00000 





00000 


-32 . 5000 





0000 





.0000 


0.00000 





00000 


-0 


24115 





79072 





00000 





00000 


-32 . 4000 





0000 





. 0000 


0.00000 





00000 


-0 


26138 





79303 





00000 





00000 


-32.3000 





0000 





. 0000 


0.00000 





00000 


-0 


28379 





79883 





00000 





00000 


-32.2000 





0000 





. 0000 


0.00000 





00000 


-0 


31030 





81135 





00000 





00000 


-32.1000 





0000 





. 0000 


0.00000 





00000 


-0 


34557 





83835 





00000 





00000 


-32.0000 





0000 





. 0000 


0.00000 





00000 


-0 


38043 





86434 





00000 





00000 


-31.9000 





0000 





. 0000 


0.00000 





00000 


-0 


41426 





88827 





00000 





00000 


-31.8000 





0000 





.0000 


0.00000 





00000 


-0 


44627 





90876 





00000 





00000 


-31.7000 





0000 





.0000 


.00000 





00000 


-0 


47560 





92430 





00000 





00000 


-31 . 6000 





0000 





.0000 


0.00000 





00000 


-0 


50149 





93347 





00000 





00000 


-31.5000 





0000 





.0000 


.00000 





00000 


-0 


52359 





93550 





00000 





00000 


-31.0000 





0000 





. 0000 


0.00000 





00000 


-0 


62995 





90101 





00000 





00000 


-30.0000 





0000 





. 0000 


0.00000 





00000 


-0 


83356 





80658 





00000 





00000 


-29 . 0000 





0000 





. 0000 


0.00000 





00000 


-1 


01850 





68353 





00000 





00000 



In the example given, the incident wave is propagating in the +x direction. The first layer of dipoles in 
the rectangular target is at x/d = —32.5, and the "surface" of the target is at x/d = —33.. We have 
evaluated the E field on the illuminated side of the target, beginning at a; = —2.5d from the first dipole 
layer, a distance 2d from the "surface" of the target, and continuing to xtf = —29d, a distance 3.5d 
into the target volume. Examination of the output file shows that the E field intensity rises slightly as we 
cross the target surface. 

28 Finale 

This User Guide is somewhat inelegant, but we hope that it will prove useful. The structure of the 
ddscat .par file is intended to be simple and suggestive so that, after reading the above notes once, 
the user may not have to refer to them again. 

Known bugs in DDSCAT wiU be posted at the DDSCAT web site, 

'http : / / www .astro. prince ton. edu /'^dr a ine/ DDSCAT . html 
and the latest version of DDSCAT will always be available at this site. There is also a wiki: 

[http : //wikidot . com/start| 
Users are encouraged to provide B. T. Draine (draine@astro . princeton . edu) with their email 
address; email notification of bug fixes and any new releases of DDSCAT, will be made known to those 
who do! 

P. J. Flatau maintains the WWW page "SCATTERLIB - Light Scattering Codes Library" with URL 
|http : / /at ol . ucsd . edu /^p flatau The SCATTERLIB Internet site is a library of light scat- 
tering codes. Emphasis is on providing source codes (mostly FORTRAN). However, other information 
related to scattering on spherical and non-spherical particles is collected: an extensive list of references 
to light scattering methods, refractive index, etc. This URL page contains section on the discrete dipole 
approximation. 

Concrete suggestions for improving DDSCAT (and this User Guide) are welcomed. To cite this User 
Guide, we suggest the following citation: 

Draine, B.T., & Flatau, RJ. 2008, "User Guide for the Discrete Dipole Approximation Code 

DDSCAT (Version 7.0)", http://arxiv.org/abs/0809.0337 
Users may also wish to cite papers describing DDA theory and its implementation, such as Draine (1988), 
Draine & Flatau (2004), and Draine & Flatau (2008). 

Finally, the authors have one special request: We would appreciate preprints and (especially!) reprints 
of any papers which make use of DDSCAT! 
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A Understanding and Modifying ddscat.par 

In order to use DDSCAT to perform the specific calculations of interest to you, it will be necessary to 
modify the ddscat . par file. Here we list the sample ddscat . par file, followed by a discussion of 
how to modify this file as needed. Note that all numerical input data in DDSCAT is read with free-format 
READ (IDEV, *) . . . statements. Therefore you do not need to worry about the precise format in which 
integer or floating point numbers are entered on a line. The crucial thing is that lines in ddscat . par 
containing numerical data have the correct number of data entries, with any informational comments 
appearing after the numerical data on a given line. 

' ========= Parameter file for v7 . . 7 ===================' 

'**** Preliminaries ****' 

'NOTORQ' = CMT0RQ*6 (NOTORQ, DOTORQ) — either do or skip torque calculations 
'PBCGS2' = CMDS0L*6 (PBCGS2, PBCGST, PETRKP) — select solution method 
'GPFAFT' = CMDFFT*6 (GPFAFT, FFTMKL) 
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR 
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN) 
'**** Initial Memory Allocation ****' 

100 100 100 = dimensioning allowance for target generation 
'**** Target Geometry and Composition 
'RCTGLPRSM' = CSHAPE*9 shape directive 
32 24 16 = shape parameters 1-3 

1 = NCOMP = number of dielectric materials 
'diel.tab' = file with refractive index 1 

Error Tolerance ****' 
l.OOe-5 = TOL = MAX ALLOWED (NORM OF | G>=AC| E>-ACA | X> ) / (NORM OF AC|E>) 
'**** Interaction cutoff parameter for PBC calculations ****' 
5.00e-3 = GAMMA (le-2 is normal, 3e-3 for greater accuracy) 
'**** Angular resolution for calculation of <cos>, etc. ****' 
0.5 = ETASCA (number of angles is proportional to [ ( 3+x) /ETASCA] ' 2 ) 

Wavelengths (micron) 
6.283185 6.283185 1 'LIN' = wavelengths ( first , last , how many , how=LIN, INV, LOG) 
'**** Effective Radii (micron) **** ' 

2 2 1 'LIN' = aeff ( first , last , how many , how=LIN, INV, LOG) 
'**** Define Incident Polarizations ****' 

(0,0) (l.,0.) (0.,0.) = Polarization state eOl (k along x axis) 

2 = lORTH (=1 to do only pol . state eOl; =2 to also do orth. pol . state) 

Specify which output files to write ****' 
1 = IWRKSC (=0 to suppress, =1 to write ".sea" file for each target orient. 

1 = IWRPOL (=0 to suppress, =1 to write ".pol" file for each (BETA, THETA) 
'**** Prescribe Target Rotations ****' 

0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around al) 

0. 90. 3 = THETMI, THETMX, NTHETA (theta=angle between al and k) 

0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of al around k) 

'**** Specify first IWAV, IRAD, lORI (normally 0) ****' 

= first IWAV, first IRAD, first lORI (0 to begin fresh) 

'**** Select Elements of S_ij Matrix to Print ****' 

6 = NSMELTS = number of elements of S_ij to print (not more than 9) 
11 12 21 22 31 41 = indices ij of elements to print 

Specify Scattered Directions ****' 
'LFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame) 

2 = NPLANES = number of scattering planes 

0. 0. 180. 10 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 1 
90. 0. 180. 10 = phi, thetan_min, thetan_max, dtheta (in deg) for plane 2 
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Lines Comments 
1-2 comment lines 

3 NOTORQ if torque calculation is not required; 
DOTORQ if torque calculation is required. 

4 PBCGS2 is recommended; other options ai-e PBCGST and PETRKP (see SjTT)- 

5 GPFAFT is supplied as default, but FFTMKL is recommended if DDSCAT has been compiled with 
the Intel® Math Kernal Library (see SC5l[T2b. 

6 GKDLDR is recommended as the prescription for polarizabilities (see ifT3] ) 

7 NOTBIN for no unformatted binary output. 

ORIBIN for unformatted binary dump of orientational averages only; 
ALLBIN for full unformatted binary dump ( ^10.21 ); 

8 comment line 

9 initial memory allocation NX,NY,NZ. These must be large enough to accomodate the target that 
will be generated. 

10 comment line 

11 specify choice of target shape (see SjT9]for description of options RCTGLPRSM, ELLIPSOID, 
TETRAHDRN, ...) 

12 shape parameters SHPARl, SHPAR2, SHPAR3, ... (see 

13 number of different dielectric constant tables (fT4l). 

14 name(s) of dielectric constant table(s) (one per line). 

15 comment line 

16 TOL = error tolerance h: maximum allowed value of \A^E - A^AP\/\A^E\ [see eq.®]. 

17 comment line 

18 GAMMA= interaction cutoff parameter 7 (see Draine & Flatau 2009) 
the value of 7 does not affect calculations for isolated targets 

19 comment line 

20 ETASCA - parameter ?/ controlling angular averages ( ^23^ . 

21 comment line 

22 A - first, last, how many, how chosen. 

23 comment line 

24 tteff - first, last, how many, how chosen. 

25 comment line 

26 specify x,y,z components of (complex) incident polarization eoi (^22l) 

27 lORTH = 2 to do both polarization states (normal); 
lORTH = 1 to do only one incident polarization. 

28 comment line 

29 IWRKSC = to suppress writing of ".sea" files; 
IWRKSC = 1 to enable writing of ".sea" files. 

30 IWRPOL = to suppress writing of ".pol" files. 
IWRPOL = 1 to write ".pol" files. 

3 1 comment line 

32 /? (see ifTTb - first, last, how many . 

33 Q (see 23 - first, last, how many. 

34 4> (see fJT7]l - first, last, how many. 

35 comment line 

36 IWAVO IRADO lORlO - Starting values of integers IWAV IRAD I ORI (normally 0). 

37 comment line 

38 Ns = number of scattering matrix elements (must be < 9) 

39 indices ij of Ns elements of the scattering matrix Sij 

40 comment line 

41 specify whether scattered directions are to be specified by CMDFRM= ' LFRAME ' or ' TFRAME ' . 

42 NP LANES = number of scattering planes to follow 

43 (j)s for first scattering plane, 9s,min, Os,max, how many 6s values; 
44,... (j)s for 2nd,... scattering plane, ... 
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B wxxxr3;_y . avg Files 

The file wOOOrOOOori. avg contains the resuhs for the first wavelength (wOOO) and first target ra- 
dius (rOO) averaged over orientations ( . avg). The wOOOrOOOori. avg file generated by the sample 
calculation should look like the following: 

DDSCAT DDSCAT 7.0.7 [08.08.29] 

TARGET Rectangular prism; NX,NY,NZ= 32 24 16 

GPFAFT method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

12288 = NATO = number of dipoles 

0.06985579 - d/aeff for this target [d=dipole spacing] 

0.139712 = d (physical units) 
AEFF= 2.000000 = effective radius (physical units) 

WAVE= 6.283185 = wavelength (physical units) 

K*AEFF= 2.000000 = 2*pi*aeff /lambda 

n= ( 1.3300 , 0.0100), eps . = ( 1.7688 , 0.0266) | m | l<d= 0.1858 for subs. 1 

TOL= l.OOOE-05 = error tolerance for CCG method 
( 1.00000 0.00000 0.00000) = target axis Al in Target Frame 
( 0.00000 1.00000 0.00000) = target axis A2 in Target Frame 

NAVG= 603 = (theta,phi) values used in comp. of Qsca,g 
( 0.13971 0.00000 0.00000) = k vector (latt. units) in Lab Frame 
( 0.00000, 0.00000) ( 1.00000, 0.00000) ( 0.00000, O.00000)=inc.pol.vec. 1 in LF 
( 0.00000, 0.00000) ( 0.00000, 0.00000) ( 1.00000, . 00000) =inc.pol.vec. 2 in LF 

0.000 0.000 = beta_min, beta_max ; NBETA = 1 

0.000 90.000 = theta_min, theta_max; NTHETA= 3 

0.000 0.000 = phi_min, phi_max ; NPHI = 1 



0.5000 = ETASCA = param 
Results averaged over 
and 

Qext 
8 . 6320E-01 
5 . 9773E-01 
7 . 3047E-01 
2 . 6547E-01 
Qsca*g ( 1 ) 
5 . 1512E-01 
3 . 7677E-01 
4 . 4595E-01 



J0=1 
J0=2 
mean 
Qpol 

J0=1 
J0=2 
mean 



controlling # of scatt . dirs used to calculate <cos> etc. 
3 target orientations 
2 incident polarizations 



Qabs 

8.1923E-02 7 

6.5842E-02 S 

7.3882E-02 f 

Qsca*g (2 ) 
2 . 0999E-02 
3 . 3543E-02 
2 . 7271E-02 



Qsca 
. 8128E- 
. 3189E- 
. 5659E- 



01 
01 
01 



g ( 1 ) =<cos> 
6 . 5932E-01 
7 . 0837E-01 
6 . 7919E-01 



<cos " 2> 
5 . 5093E-01 
5 . 7263E-01 
5 . 7263E-01 



. 6216E-03 
. 6332E-03 
. 6274E-03 
dQpha= 



Qpha 

9 . 5911E-01 
8 . 9157E-01 
9 . 2534E-01 
6.7534E-02 



Qsca*g ( 3 ) 
2 . 0275E-08 
-4 . 5799E-09 
7 . 8476E-09 



iter 
3 
3 



mxiter 
13 
13 



Nsca 
603 
603 



Mueller matrix elements for selected scattering directions in Lab Frame 



theta 


phi 




Pol . 


S_ll 




3_12 




S_21 




S_22 


S_31 




S_41 







00 





00 





11122 


3 


9825E+00 


4 


. 4295E- 


01 


4 


430E- 


01 


3 


982E+00 


-9 


951E- 


09 


-1 


167E- 


08 


10 


00 





00 





09797 


3 


7743E+00 


3 


. 6976E- 


01 


3 


698E- 


01 


3 


774E+00 


1 


511E- 


09 


-8 


158E- 


09 


20 


00 





00 





06514 


3 


2121E+00 


2 


. 0924E- 


01 


2 


092E- 


01 


3 


212E+00 


8 


423E- 


09 


-9 


677E- 


09 


30 


00 





00 





01115 


2 


4761E+00 


2 


. 7618E- 


02 


2 


762E- 


02 


2 


476E+00 


6 


078E- 


08 


1 


683E- 


08 


40 


00 





00 





06519 


1 


7463E+00 


-1 


. 1383E- 


01 


'1 


138E- 


01 


1 


746E+00 


-1 


787E- 


08 


2 


854E- 


08 


50 


00 





00 





16289 


1 


1386E+00 


-1 


. 8546E- 


01 


-1 


855E- 


01 


1 


139E+00 


-9 


912E- 


09 


-4 


191E- 


08 


60 


00 





00 





27512 


6 


9235E-01 


-1 


. 9048E- 


01 


'1 


905E- 


01 


6 


923E-01 


-1 


481E- 


08 


-1 


639E- 


08 


70 


00 





00 





38377 


3 


9361E-01 


-1 


. 5106E- 


01 


-1 


511E- 


01 


3 


936E-01 


2 


248E- 


08 


-6 


877E- 


09 


80 


00 





00 





45466 


2 


0690E-01 


-9 


. 4071E- 


02 


-9 


407E- 


02 


2 


069E-01 


1 


167E- 


09 


-2 


473E- 


08 


90 


00 





00 





43502 


9 


7450E~02 


-4 


.2393E- 


02 


-4 


239E- 


02 


9 


745E-02 


-3 


071E- 


09 


4 


950E- 


09 


100 


00 





00 





24195 


4 


0410E-02 


-9 


. 7771E- 


03 


-9 


777E- 


03 


4 


041E-02 


-5 


077E- 


09 


-2 


816E- 


08 


110 


00 





00 





09530 


2 


0545E-02 


1 


. 9580E- 


03 


1 


958E- 


03 


2 


054E-02 


1 


702E- 


09 


-1 


171E- 


08 


120 


00 





00 





03721 


2 


7014E-02 


1 


. 0051E- 


03 


1 


005E- 


03 


2 


701E-02 


1 


912E- 


10 


7 


689E- 


10 


130 


00 





00 





03732 


4 


8444E-02 


-1 


. 8079E- 


03 


-1 


808E- 


03 


4 


844E-02 


-1 


631E- 


09 


-3 


689E- 


09 


140 


00 





00 





00713 


7 


2647E-02 


-5 


. 1804E- 


04 


-5 


180E- 


04 


7 


265E-02 


2 


516E- 


09 


1 


104E- 


09 


150 


00 





00 





04843 


9 


0620E-02 


4 


. 3883E- 


03 


4 


388E- 


03 


9 


062E-02 


-8 


039E- 


10 


2 


467E- 


09 


160 


00 





00 





09683 


9 


9496E-02 
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10 
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2 
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C ^fjxxxryyykzzz .sea Files 



The wOOOrOOOkOOO . sea file contains the resuhs for the first wavelength (wOOO), first target radius 
(rOO), and first orientation (kOOO). The wOOOrOOOkOOO . sea file created by the sample calculation 
should look like the following: 

DDSCAT DDSCAT 7.0.7 [08.08.29] 

TARGET Rectangular prism; NX,NY,NZ= 32 24 16 

GPFAFT method of solution 

GKDLDR prescription for polarizabilies 

RCTGLPRSM shape 

12288 = NATO = number of dipoles 
0.06985579 = d/aeff for this target [d=dipole spacing] 
0.139712 = d (physical units) 
physical extent of target volume in Target Frame 



-4.470771 0. 
-1.676539 1. 
-1.117693 1. 
AEFF= 2.000000 
WAVE= 6.283185 
K*AEFF= 2.000000 
n= ( 1.3300 , 0.0100) 



000000 = xmin,xmax (physical units) 

676539 = ymin,ymax (physical units) 

117693 = zmin, zmax (physical units) 

= effective radius (physical units) 
= wavelength (physical units) 
= 2*pi*aeff /lambda 

eps.= ( 1.7688 , 0.0266) lmll<d= 



0.1858 for subs. 1 



TOL= l.OOOE-05 
1.00000 0.00000 



.00000 
NAVG= 
.13971 
00000, 
00000, 
00000 
00000 
13971 
00000, 
00000, 



error tolerance for CCG method 
0.00000) = target axis Al in Target Frame 
1.00000 0.00000) = target axis A2 in Target Frame 
603 = (theta,phi) values used in comp. of Qsca,g 



0.00000 0.00000) = 
0.00000) ( 1.00000, 
0.00000) ( 0.00000, 
0.00000 0.00000) = 
1.00000 0.00000) = 
0.00000 0.00000) = 
0.00000) ( 1.00000, 



k vector (latt. units) in TF 
0.00000) ( 0.00000, . 00000) =inc .pel .vec . 
0.00000) ( 1.00000, O.00000)=inc.pol.vec. 
target axis Al in Lab Frame 
target axis A2 in Lab Frame 
k vector (latt. units) in Lab Frame 
0.00000) ( 0.00000, . 00000) =inc .pol .vec . 

0. 00000 )=inc.pol. vec. 



1 in TF 

2 in TF 



1 in LF 

2 in LF 



0.00000) ( 0.00000, 0.00000) ( 1.00000, 
BETA = 0.000 = rotation of target around Al 
THETA= 0.000 = angle between Al and k 

PHI = 0.000 = rotation of Al around k 
0.5000 = ETASCA = param. controlling # of scatt . dirs used to calculate <cos> etc. 



J0=1 
J0=2 
mean 
Qpol 

J0=1 
J0=2 
mean 



Qpha 
9.7717E- 
9. 1277E- 
9.4497E- 
6.4397E- 



Qext Qabs Qsca g(l)=<cos> <cos"2> Qbk 

9.0997E-01 8.6105E-02 8.2387E-01 6.4090E-01 5.8427E-01 2.5108E-02 

6.9891E-01 7.0519E-02 6.2839E-01 6.5894E-01 6.0509E-01 1.9097E-02 

8.0444E-01 7.8312E-02 7.2613E-01 6.4871E-01 5.9328E-01 2.2103E-02 

2.1106E-01 dQpha= 

Qsca*g(l) Qsca*g(2) Qsca*g(3) iter mxiter Nsca 

5.2802E-01 -1.2657E-08 -1.6907E-08 3 13 603 

4.1407E-01 -4.0495E-08 -1.0475E-09 3 13 603 

4.7105E-01 -2.6576E-08 -8.9774E-09 

Mueller matrix elements for selected scattering directions in Lab Frame 
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3 


. 7772E- 


02 


3 


. 777E- 


02 


2 


. 778E-01 
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30 


00 


90 


00 





20414 


3 


2267E+00 


-6 


5869E- 


01 


-6 


587E- 


01 


3 


227E+00 


-2 


834E- 


08 


-9 


660E- 


10 


40 


00 


90 


00 





29034 


2 


5486E+00 


-7 


3995E- 


01 


-7 


400E- 


01 


2 


549E+00 


-1 


776E- 


08 


1 


966E- 


08 


50 


00 


90 


00 





40271 


1 


8295E+00 


-7 


3678E- 


01 


-7 


368E- 


01 


1 


830E+00 


-6 


937E- 


09 


4 


586E- 


09 


60 


00 


90 


00 





53656 


1 


1717E+00 


-6 


2868E- 


01 


-6 


287E- 


01 


1 


172E+00 


-1 


026E- 


09 


1 


203E- 


08 


70 


00 


90 


00 





67801 


6 


5188E-01 


-4 


4198E- 


01 


-4 


420E- 


01 


6 


519E-01 


-2 


419E- 


09 


6 


202E- 


09 


80 


00 


90 


00 





79934 


2 


9924E-01 


-2 


3919E- 


01 


-2 


392E- 


01 


2 


992E-01 


7 


830E- 


10 


7 


416E- 


09 


90 


00 


90 


00 





85269 


9 


9116E-02 


-8 


4516E- 


02 


-8 


452E- 


02 


9 


912E-02 


5 


458E- 


10 


1 


582E- 


09 


100 


00 


90 


00 





60400 


1 


5341E-02 


-9 


2659E- 


03 


-9 


266E- 


03 


1 


534E-02 


-4 


916E- 


10 


9 


156E- 


10 


110 


00 


90 


00 





19961 


1 


0545E-02 


-2 


1049E- 


03 


-2 


105E- 


03 


1 


055E-02 


1 


321E- 


10 


-6 


863E- 


10 


120 


00 


90 


00 





51133 


5 


3315E-02 


-2 


7261E- 


02 


-2 


726E- 


02 


5 


332E-02 


3 


564E- 


10 


-1 


670E- 


09 


130 


00 


90 


00 





44323 


1 


1701E-01 


-5 


1864E- 


02 


-5 


186E- 


02 


1 


170E-01 


9 


436E- 


10 


-6 


678E- 


10 


140 


00 


90 


00 





34103 


1 


7989E-01 


-6 


1350E- 


02 


-6 


135E- 


02 


1 


799E-01 


-2 


355E- 


10 


1 


092E- 


10 


150 


00 


90 


00 





25177 


2 


2835E-01 


-5 


7491E- 


02 


-5 


749E- 


02 


2 


283E-01 


2 


310E- 


09 


4 


675E- 


10 


160 


00 


90 


00 





18698 


2 


5853E-01 


-4 


8340E- 


02 


-4 


834E- 


02 


2 


585E-01 


4 


504E- 


09 


5 


089E- 


10 


170 


00 


90 


00 





14861 


2 


7346E-01 


-4 


0639E- 


02 


-4 


064E- 


02 


2 


735E-01 


3 


817E- 


09 


3 


120E- 


10 


180 


00 


90 


00 





13599 


2 


7775E-01 


~3 


7772E- 


02 


-3 


777E- 


02 


2 


778E-01 


3 


076E- 


09 


-1 


357E- 


10 



WXXXRYYYKZZZ .POLN HLES 



D ^fjxxxryyykzzz .poln Files 

The wOOOrOOOkOOO.poll file contains the polarization solution for the first wavelength ( w ), first 
target radius (r 0), first orientation (kO 0), and first incident polarization (pol 1). In order to limit the 
size of this file, it has been written as an unformatted file. This preserves full machine precision for the 
data, is quite compact, and can be read efficiently, but unfortunately the file is not fully portable because 
different computer architectures (e.g., Linux vs. MS Windows) have adopted different standards for 
storage of "unformatted" data. However, anticipating that many users will be computing within a single 
architecture, the distribution version of DDSCAT uses this format. 

Additional warning: even on a single architecture, users should be alert to the possibility that different 
compilers may follow different conventions for reading/writing unformatted files. 

The file contains the following information: 

• The location of each dipole in the target frame. 

• {kx, ky, kz)d, where k is the incident k vector. 

• (£^02; , Eoy , Eoz), the complex polarization vector of the incident wave. 

• a~^d^, the inverse of the symmmetric complex polarizability tensor for each of the dipoles in the 
target. 

• {Px, Py, Pz), the complex polarization vector for each of the dipoles. 

The interested user should consult the routine writepol . f, or the to see how this information has 
been organized in the unformatted file. The user can, of course, simply use the routine readpol . f to 
read this file and load the information into memory. 

Note that a subroutine readpol . f 90 is already provided that reads the stored file, and the program 
DD field. f90 (which uses readpol . f 90) is provided to calculate E and B at user-selected postions 
-seefgT] 
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Qabs, 8 
Qcxt, 8 
Qsca.^ 8 
Qpha, 8 

Sj - scattering amplitude matrix, 59 

Stj - 4x4 Mueller scattering matrix, 59, 60, 62 

Sj - scattering amplitude matrix, 60 

$ - target orientation angle, 20, 30 

<i>DF, 40, 64 

8 - target orientation angle, 20, 30 
Gdf, 40, 64 

P - target orientation angle, 20, 30 
PuF, 40, 64 

7] - parameter for choosing scattering angles, 59 
7] - parameter for selection of scattering angles, 57 
ai, 3.2 target vectors, 20 
(ps - scattering angle, 52, 57 
9s - scattering angle, 52, 57 
Ocff , 7, 20 

fij - amplitude scattering matrix, 59 
fij - scattering amplitude matrix, 60 
Qr - radiative torque efficiency vector, 28 
Qpr - radiative force efficiency vector, 28 

absorption efficiency factor Qabs, 8 
averages over scattering angles, 56 

CALLTARGET, 41 

CALLTARGET and PBC, 42 

CMDFFT - specifying FFT method, 70 

CMDFRM - specifying scattering directions, 52, 70 

compiling and linking, 14 

conjugate gradient algorithm, 24 

Cygwin, 14, 18 

DDfield, 20, 65 

ddscat.log.OOO - output file, 23 

ddscat.out, 23 

ddscat.par 

ALLBIN, 24 

CSHAPE, 32 

diel.tab, 19 

DOTORQ, 28 

ETASCA, 28 

GAMMA, 19 

GKDLDR, 19, 26 

GPFAFT, 19 

lORTH, 59, 63 

IWRKSC, 23 

LATTDR, 19, 26 

NOTBIN, 19 

ORIBIN, 24 



PBCGST, 24 
PETRKP, 24 
SHPARi,SHPAR2,..., 32 
TOL, 24 

ddscat.par - BETAMI,BETAMX,NBETA, 31 
ddscat.par - PHIMIN,PHIMAX,NPHI, 31 
ddscat.par - THETMI,THETMX,NTHETA, 3 1 
ddscat.par parameter file, 18 
diel.tab - see ddscat.par, 19 
Dielectric Frame (DF), 34, 40, 64 
dielectric function of target material, 26 
dielectric medium, 1 1 
dielectric tensor, 64 

differential scattering cross section dCsca/dil, 62 
DOTORQ, 28 

effective radius Ocff , 7, 20 

Electric field in or near the target, 65 

ELLIPSOID, 37 

error tolerance, 24 

ESELF 68 

ETASCA, 28 

extinction efficiency factor Qcxt, 8 

FFT algorithm, 25 
FFTW, 25 
Fortran compiler 

optimization, 15 

GAMMA - see ddscat.par, 19 
GKDLDR, 26 

GKDLDR - see ddscat.par, 19 
GPFA, 25 
GPFAFT, 47 

GPFAFT - see ddscat.par, 19 
GPFAPACK, 68 

IDVERR, 16 

IDVOUT, 16 

infinite cylinder, 43 

infinite hexagonal column, 47 

lORTH, 63 

lORTH parameter, 59 

iterative algorithm, 24 

IWRKSC - see ddscat.par, 23 

IWRPOL, 65 

Lab Frame (LF), 29 
LATTDR, 26 

LATTDR - see ddscat.par, 19 
LFRAME, 52 
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Magnetic field in or near the target, 65 

memory requirements, 28 

MKL 

Intel® MKL 
DFTI, 15 

MPI - Message Passing Interface, 16 

code execution, 22 
MS-Windows and other abominations, 17 
mtable - output file, 23 

Mueller matrix for infinite targets, periodic in 1-d, 63 
Mueller matrix for infinite targets, periodic in 2-d, 63 
Mueller matrix for scattering, 59, 60, 62 
MXNX,MXNY,MXNZ, 28 

NCOMP 

ddscat.par, 19 
NOTBIN - see ddscat.par, 19 

orientational averaging, 3 1 

nonrandomly-oriented targets, 32 

randomly-oriented targets, 32 
orientational sampling in (3, 8, and $,31 

PBC = periodic boundary conditions, 54 
PBCGST algorithm, 24 

periodic boundary conditions, 36, 42, AA-Al, 49-5 1 

PETRKP algorithm, 24 

phase lag efficiency factor Qpha, 8 

PIM package, 24 

polarizabilities 

GKDLDR, 26 

LATTDR, 26 
polarization - elliptical, 59 
polarization of incident radiation, 55 
polarization of scattered radiation, 61, 62 
precision: single vs. double, 15 
PYD, 42, 46, 53, 54 
PZD, 42, 46, 54 

qtable - output file, 23 
qtable2 - output file, 23 
qtable2 file, 23 

radiative force efficiency vector Qp,., 28 
radiative torque efficiency vector Qr, 28 
READPOL, 65 

refractive index of target material, 26 
relative dielectric function, 1 1 
relative refractive index, 1 1 

scattering - angular averages, 56 
scattering angles Os, (j>s, 57 
scattering by tilted cube, 56 
scattering directions, 52 
scattering efficiency factor Qsca, 8 



SCATTERLIB, 67 
size parameter x = fcocff , 7 
source code (downloading of), 13 
Stokes vector (/, Q, U, V), 60 

Target Frame, 64 
Target Frame (TF), 30 
target generation, 32 

target generation: infinite periodic targets, 42 
target orientation, 29, 30 
target routines: modifying, 41 
target shape options 

AN1_ELL_2, 34 

ANI_ELL_3, 34 

ANIELLIPS, 34 

ANIFRMFIL, 33 

ANIRCTNGL, 35 

BISLINPBC, 42 

CONELLIPS, 36 

CYLINDERl, 36 

CYLNDRCAP 36 

CYLNDRPBC, 42 

DSKBLYPBC, 44 

DSKRCTNGL, 36 

DSKRCTPBC, 45 

DW1996TAR, 37 

ELLIPSO^, 37 

ELLlPSO_3, 37 

ELLIPSOID, 37 

FROM_FILE, 33 

HEX J-RISM, 38 

HEXGONPBC, 46 

LAYRDSLAB, 38 

LYRSLBPBC, 47 

MLTBLOCKS, 38 

RCTGL_PBC, 49 

RCTGLPRSM, 38 

RECRECPBC, 50 

SPH_ANI_N, 40 

SPHERES.N, 39 

SPHRN.PBC, 51 

SPHROID_2, 40 

TETRAHDRN, 41 

TRILYRPBC, 52 

TRNGLPRSM, 41 

UNIAXICYL, 41 
target shape optionsz 

RCTGLBLK3, 39 
target.out - output file, 42 
TFRAME, 52 

wOOOrOOOkOOO.poln - output file, 23 
wOOOrOOOkOOO.sca- output file, 23 
wOOOrOOOori.avg- output file, 23 



INDEX 



wOOOrOOOori.avg file - output file, 
winzip, 14 



